A dynamic discrete dislocation plasticity method for the simulation of plastic relaxation under shock loading

In this article, it is demonstrated that current methods of modelling plasticity as the collective motion of discrete dislocations, such as two-dimensional discrete dislocation plasticity (DDP), are unsuitable for the simulation of very high strain rate processes (106 s−1 or more) such as plastic relaxation during shock loading. Current DDP models treat dislocations quasi-statically, ignoring the time-dependent nature of the elastic fields of dislocations. It is shown that this assumption introduces unphysical artefacts into the system when simulating plasticity resulting from shock loading. This deficiency can be overcome only by formulating a fully time-dependent elastodynamic description of the elastic fields of discrete dislocations. Building on the work of Markenscoff & Clifton, the fundamental time-dependent solutions for the injection and non-uniform motion of straight edge dislocations are presented. The numerical implementation of these solutions for a single moving dislocation and for two annihilating dislocations in an infinite plane are presented. The application of these solutions in a two-dimensional model of time-dependent plasticity during shock loading is outlined here and will be presented in detail elsewhere.


Introduction (a) Plasticity and dislocation dynamics
Plastic deformation in crystalline materials occurs through the motion of defects. Point defects may be involved as in diffusional creep, and interfaces may be involved as in twinning and stressinduced martensitic transformations. Plastic deformation through the generation and motion of dislocations alone shall be the sole concern of this article.
The aim of the technique known as discrete dislocation dynamics (DDD) is to simulate plasticity as the result of the collective motion of individual dislocations within an elastic continuum. That is, dislocations are modelled as Volterra singularities in the continuum, and plasticity arises as the result of their generation and movement. Long-range interactions between dislocations are mediated by their elastic fields. Short-range interactions, such as dislocation reactions, pinning by obstacles and annihilations, may be included through the introduction of rules that apply when dislocations come within a certain distance of one another. The underlying theoretical framework is provided by the theory of elasticity, and both isotropic and anisotropic elastic formulations of DDD exist [1]. The Burgers vectors of the dislocations and their slip planes are determined by the crystal structure.
Furthermore, it is not enough to describe the elastic state arising from the dislocation microstructure [2]. The second essential ingredient in a treatment of DDD is the force that makes a dislocation segment move, and how this force relates to the motion of the dislocation segment itself. A segment of dislocation line will experience a force whenever its movement lowers the potential energy of the system. The potential energy comprises the total elastic energy and the potential energy associated with the loads applied to the surface of the body. The force is therefore conservative and it is independent of the description of the distortion of the medium provided by elasticity. The force is given by the usual Peach-Koehler expression [3], where σ is the stress tensor acting on the segment, ξ is the unit tangent vector of the dislocation segment, B the Burgers vector and '∧' denotes the vector product. Mura [4] showed that the Peach-Koehler expression remains valid even when the dislocation is moving, and there is no equivalent to a Lorentz force acting on a moving dislocation. When a dislocation moves through a crystal lattice, it excites atomic vibrations, and phonons are scattered by the dislocation [5,6]; at higher speeds, there may also be significant electronic excitations [6]. The friction that these interactions creates leads to severe over-damping of the dislocation motion, with the result that a dislocation quickly reaches a terminal glide velocity v [5]. When dislocations move at modest speeds in relation to the transverse wave speed, it is observed that the damping force is proportional to the glissile component of the Peach-Koehler force f PK acting on the corresponding dislocation segment [5,7,8], where d is the viscous drag coefficient of the dislocation. The physical origin and value of d are topics of ongoing research. The authors are aware of the limitations of equation (1.2), something that will be addressed in a future paper.

(b) Discrete dislocation plasticity
Discrete dislocation plasticity (DDP) is a two-dimensional formulation of DDD, describing plasticity as the collective quasi-static motion of dislocations modelled as line singularities in a linear elastic solid. This necessarily limits the scope of the method to infinite straight edge dislocations so that plane strain conditions apply. The DDP method proposed by Van der Giessen & Needleman [9] enabled interesting questions associated with dislocation-mediated plasticity and failure to be investigated, such as size effects in plastic flow [10,11], geometrical effects [12], fracture mechanics [13][14][15], crack growth [16], fatigue [17] and more. Thus, despite the obvious shortcomings (e.g. lack of cross-slip and junction formation) of a two-dimensional model of plasticity, DDP has provided useful insight into a wide range of problems. The quasi-static approximation. In the DDP methodology, as in most dislocation dynamics methods, there is an implicit assumption that the elastic field of a moving dislocation is propagated instantaneously; in DDD, the term dynamics is used to imply time evolution of the dislocation structure, rather than elastodynamics. This assumption enables the stress field at an instant in time to be evaluated by considering the static elastic fields of the dislocations at their current positions. This is reasonable as long as the speeds of dislocations are a small fraction of the elastic transverse speed of sound, which is the slowest speed at which information can be propagated in an elastic medium.
It has been speculated that at high enough strain rates, the quasi-static approximation breaks down (cf. [5,18,19]). A simple estimate shows that this is a reasonable expectation. Consider Orowan's equation [20],˙ = Bρ mv , (1.3) where˙ is the macroscopic strain rate, B the magnitude of the Burgers vector, ρ m the mobile dislocation density andv the average velocity of mobile dislocations. With ρ m ≈ 10 13 m −2 , a strain rate of 10 6 s −1 and B ≈ 0.2 nm, Orowan's equation predicts an average dislocation velocity of the order of 10 3 m s −1 , which is the same order of magnitude as the transverse speed of sound in many metals. However, it has often been argued that the limiting speed of edge dislocations is the transverse speed of sound. As Frank [21], Eshelby [22] and Weertman [23], amongst others, have pointed out, the elastic fields of uniformly moving edge dislocations diverge at that speed, making the associated elastic energy infinite. Thus, the transverse speed of sound cannot be reached by a dislocation according to the theory of linear, local, first-order elasticity. Furthermore, as a dislocation approaches the transverse speed of sound, its elastic field contracts in the direction of the motion. This is comparable with the Fitzgerald contraction of the electric field of a moving charge in the direction of its motion as its speed approaches that of light [24]. This contraction of the elastic field is usually referred to as a dislocation inertia, or a dynamic or relativistic effect.
The need for a dynamic formulation. Orowan's equation suggests that these relativistic effects may be expected at strain rates of the order of 10 6 s −1 . At those strain rates, the shapes of the elastic fields of moving dislocations may differ from their static counterparts and, thus, the plastic response of the material can be different. However, most DDP simulations consider situations where the strain rate does not exceed 10 1 -10 4 s −1 [25,26]. In fact, higher strain rates are rarely achieved outside shock loading, where strain rates of the order of 10 5 -10 8 s −1 are typical [27]. Hence, the quasi-static approximation is justified in most DDP applications.
A number of DDD studies at high strain rates exist. For instance, Roos et al. [25,26] adapted the usual DDP methodology to include the effect of relativistic contractions on the plastic response of materials subjected to high strain shearing loads.
Hirth et al. [28] introduced a scheme to treat dislocation dynamics beyond the quasi-static approximation in the sense that the radiation of energy by an accelerating dislocation was considered to confer on the dislocation an effective mass and an inertial force. Zbib et al. [29] incorporated the relativistic mass of Hirth et al. [28] in their method of three-dimensional DDD, which was further expanded by Zbib & Diaz de la Rubia [30]. This method has been used in a number of studies of high strain rate and shock loading plasticity, such as those by Shehadeh et al. [31] and Shehadeh [32], who have investigated the effects of strain rate, nucleation mechanisms and shock front rise times over different crystalline structures.
This work focuses on the time dependence of the evolving elastic fields caused by the generation and non-uniform motion of straight edge dislocations at speeds up to the transverse speed of sound. These fields make it possible to incorporate inertial and retardation effects on The simulated system: a two-dimensional rectangular block shocked with a high-pressure load on one end, with the other modelled as a reflective boundary. This condition is unrealistic, but serves the purpose of the simulation in making the system well defined; the simulation is stopped when the front arrives at the reflective boundary, and that is known in advance: for aluminium and a sample length of 40 µm and width 20 µm, the normal front propagates at c n = 6273 m s −1 , taking t = 6.38 ns to reach the far end.
the Peach-Koehler forces in simulations of plastic relaxation under shock loading. The term 'dynamic discrete dislocation plasticity' (D3P) has been coined for this new way of simulating dislocation dynamics. The paper focuses on the practical implementation of D3P in a numerical simulation as well as the derivations of the time-dependent elastic fields. The structure of the paper is as follows. In §2 it is shown that the quasi-static formulation of the response of dislocations to shock loading leads to the unphysical result of dislocation sources being activated ahead of the shock front. This provides the primary motivation for the development of D3P. In §3, building on the work of Markenscoff & Clifton [33], the elastodynamic fields for an injected, non-uniformly moving straight edge dislocation are derived. Section 4 expands on the implementation of these fields, and §5 shows some numerical applications of the theory.

Spurious dislocation generation ahead of the shock front
Typically, the experiments used to study the behaviour of materials at high strain rates involve a high-velocity impact between a pair of thin plates. This impact produces a state of uniaxial strain at the loading surface, simplifying the determination of the shocked state through consideration of the one-dimensional conservation equations. The dimensions of these plates are chosen to prevent unloading events originating at the plate periphery from being communicated to the central diagnosed portion of the target.
In the present work, the suitably termed 'parallel-plate impact experiment' has been simplified by considering the target plate as a two-dimensional rectangular block, and substituting the impact with a high-pressure boundary condition on the loading edge. As shown in figure 1, the 40 × 20 µm two-dimensional block was shocked with a high-pressure dynamic load under plane strain conditions, the resulting wavefront propagating through the material at the longitudinal speed of sound. It is acknowledged that these dimensions do not represent the usual aspect ratio of the specimens, but it is believed that this will not affect the results presented here.
In a first attempt to simulate the experiment, the traditional quasi-static DDP methodology proposed by Van der Giessen & Needleman [9] was used. The method first considers dislocations as mobile defects in an infinite medium. As discussed in §3f, the finite size of the simulated system is treated by exploiting the linear superposition principle [9,34] to include fields that ensure the surface of the body is free of tractions.
To simulate the dynamic loading of the two-dimensional block, the left surface was loaded with a constant uniaxial stress at t = 0, exciting an elastodynamic wavefront propagating through the solid; the upper and lower surfaces were defined as traction-free surfaces, whereas over the right surface, a reflective boundary condition was applied. A commercial finite-element package, ABAQUS, was coupled to the DDP code via Python scripting, and an explicit solver  was used for the solution of the associated elastodynamic, finite-size problem. As in usual DDP simulations, the fields of dislocations were assumed to be elastostatic and therefore to propagate instantaneously.
The elastic parameters of aluminium were used 1 in the simulation. Slip planes were oriented at ±45 • and 90 • to the direction of impact, and spaced by 100 Burgers vectors. A random population of sources and obstacles was assumed in the slip planes with a density of ρ s = 100 µm −2 . Motion of dislocations was assumed to be overdamped, with a viscous drag coefficient d = 10 −5 Pa s. A mobility law as in equation (1.2) was assumed, and the speed of dislocations was capped at the shear speed of sound. A forward Euler integration scheme was used, with a dynamic time step that limited dislocation motion to 1 nm per time step. In this first approach, heterogeneous nucleation alone was considered; hence, only low-intensity sources were included, with a strength of 100 ± 10 MPa taken from a Gaussian distribution, and the obstacle strength was set at 100 MPa. Figure 2a and b shows the positions of the wavefront and dislocation configurations obtained at t = 0.9 ns and t = 2 ns, respectively. Dislocations are seen to nucleate ahead of the front as a result of the stress fields originating from the dislocations generated behind the front. Because of the quasi-static approximation, the elastic fields of dislocations are transmitted throughout the sample at the instant the dislocations are created behind the front. This is completely unphysical, and is a direct consequence of the quasi-static assumption. In reality, these stress fields would take a finite time to be propagated, and therefore it would not be possible to activate dislocation sources until the elastic front has reached them, i.e. plastic deformation cannot be propagated faster than elastic deformation. Figure 2 indicates that these incorrect consequences of the quasistatic approximation will be avoided by solving the elastodynamic equations rather than the elastostatic equations for dislocation generation, annihilation and motion in simulations of shock loading.

Dynamic discrete dislocation plasticity
In §2, it has been argued that if dislocation-mediated plasticity is to be used for the simulation of rapid dynamic events such as shock loading, the elastic fields of dislocations have to be treated dynamically as well. In a two-dimensional model, this entails the use of dynamic solutions for the injection of straight edge dislocations and for their glide at non-uniform speeds.
It will be shown that the use of dynamic solutions of the elastic fields of dislocations fundamentally alters the DDP paradigm: all elastic interactions, instantaneous in the quasi-static approximation, change to those based on a retardation principle akin to that observed for electric charges moving at a large fraction of the speed of light in electrodynamics.

(a) Elastic waves in plane strain
Any elastodynamic problem is governed by the equation of conservation of linear momentum. In isotropic elasticity, the conservation equation can be written in the usual form of the Navier-Lamé equation [35,36], where hereafter a repeated index denotes summation, Λ and μ are Lamé's first and second constants, and ρ is the density of the medium; u i is the ith component of the elastic displacement vector, a function of position and time; u i,j denotes the first partial derivative of u i with respect to It is well known that equation (3.1) can be separated into two separate wave equations by expressing the displacement vector as the sum of the gradient of a scalar potential and the curl of a vector potential [35,36]. In the (x,z)-plane under plane strain conditions, with u y = 0 and ∂/∂y( ) ≡ 0, this process results in two wave equations in the scalar potentials φ = φ(x, z, t) and ψ = ψ(x, z, t), and In terms of these scalar potentials, the components of the displacement vector become and

Injection of a non-uniformly moving edge dislocation
In this section, the elastodynamic solution for a non-uniformly moving injected edge dislocation is developed. The method and notation introduced in 1981 by Markenscoff & Clifton [33] are used. These authors obtained the solution for a pre-existing straight edge dislocation moving at a non-uniform speed. Their solution was based on the method developed by Markenscoff [37] for the non-uniform motion of a screw dislocation. Brock [38] produced an equivalent solution to that of Markenscoff & Clifton for a pre-existing straight edge dislocation moving at a non-uniform speed. A procedure similar to that of Markenscoff was used by Jokl et al. [39] to solve the injection of a static screw dislocation. Recently, Pellegrini [40] has expanded the Peierls-Nabarro equations into a fully dynamic formulation. Consider an infinite straight edge dislocation moving in the x-direction, whose line is parallel to the y-axis as depicted in figure 3. Its position with respect to the origin at a given instant in time t is defined by a piecewise continuous function l(t), called the past history function.
The solution for the injection of a non-uniformly moving edge dislocation can be thought of as the superposition of two contributions.
-The injection contribution, describing the injection of a static edge dislocation with Burgers -The mobile contribution, describing the non-uniform motion of an existing edge dislocation, the solution for which was obtained by Markenscoff & Clifton [33], Here, H(x) is the Heaviside step function. From equation (3.6), when the dislocation is injected at t = 0, slip by the Burgers vector takes place in the semi-infinite plane z = 0, x < 0. The slipped region of the plane z = 0 is extended at t > 0 by the strip 0 ≤ x ≤ l(t) through the subsequent motion of the dislocation in equation (3.7). A further boundary condition is given by which must be fulfilled in all cases. This ensures that there is no normal stress anywhere on the slip plane at any time arising from the injection and subsequent motion of the dislocation. Solution procedure. The solution procedure for both the injection and the mobile boundaryvalue problems is based on the Cagniard-de Hoop technique [41,42]. Accordingly, first define the following Laplace transform in time and the bilateral Laplace transform in x, and where s appears in the exponential as a scaling factor for convenience. Successively applying these transforms, and assuming that at t = 0 the dislocation is at rest, equations (3.2) and (3.3) are transformed into and The solutions to equations (3.11) and (3.12) must be of the form (3.14) The integration constants C(λ, s) and C (λ, s) are found by satisfying the transformed boundary conditions for each of the contributions.
The injection contribution is the only case that shall be studied here. Transforming equations (3.6) and (3.7), the following expressions for the transformed potentials are obtained: and The inversion of the transformed solution is performed by applying the Cagniard-de Hoop technique [41,42].The procedure is illustrated for the shear component of stress, σ xz , Equation (3.17) can be transformed to Using equations (3.15) and (3.16), Let the following derivation serve as an example of how the Cagniard-de Hoop technique is applied: consider, for instance, the normal contribution-i.e. the first term in equation (3.19), The inverse bilateral Laplace transform is defined as Note the presence of the scaling factor s in the integrand. Applying it to equation (3.20), The Cagniard form of this integral must now be found; i.e. equation (3.22) must be rewritten as a forward Laplace transform by making a suitable change of integration variable and, concurrently,  of the integration path. Consider the following change of integration variable: where τ ≥ 0. τ can be expressed as It follows that λ and α are, with respect to τ , The contour of integration in the λ-plane is a Bromwich contour along the imaginary λ-axis. Upon changing variable to τ , the Bromwich contour can be distorted into a hyperbolic path as shown in figure 4. Indeed, invoking Cauchy's theorem and Jordan's lemma [33], the integral in the λ-plane alongside the Bromwich contour is seen to be equivalent to the one in the same λ-plane along the hyperbola in figure 4. The latter corresponds to an integral in the τ -plane between τ = ra and τ → ∞.
Hence, the Cagniard form of the integral is obtained, where H(τ − ra) is a Heaviside function and τ − ra a retarded time.
Applying now the inverse Laplace transform in time, it can be clearly seen that because the inverse Laplace transform is applied over an expression written in the form of a forward Laplace transform, the solution is the integrand itself,  that is, The same can be done for the rest of the terms and components. The results are summarized in table 1.

(c) Asymptotic behaviour of the static injected solution
It can be checked that the injection contributions tend to the traditional quasi-static solution in the t → ∞ limit. For instance,

Since
where B = u/2 is the Burgers vector, the injection contribution is seen to converge to the static one. The same can be proved for all other components of stress and displacement.

(d) The mobile contributions
The derivation of the mobile contributions [33] is analogous to that of the injected ones, but with the added complication that, in this case, the boundary condition depends on a past history function l(t) that is, in general, unknown.
The transformed governing equations are the same as before (equations (3.13) and (3.14)), but the integration constants must now fulfil the boundary condition above, which is first written in its equivalent form as is the inverse of the past history function) and then transformed accordingly, to finally obtain the following transformed potentials: The inversion of these potentials is performed using the Cagniard-de Hoop method following the same procedure as before. The resulting expressions are given in table 2.

(e) The superimposed solution
Section 3b summarizes the derivations of the elastodynamic fields for the injection of a static and a non-uniformly moving straight edge dislocation. The global solution for the motion of an injected straight edge dislocation is then obtained by superposition of the elastic fields associated with the mobile and injected solutions. Furthermore, the present formulation can be easily accommodated to that of a pre-existing dislocation (derived by [33]): the injection contributions would be substituted for the corresponding elastostatic solutions [18,43].

(f) The solution scheme: image fields as a linear superposition
The expressions of the elastic fields presented in the previous section are valid only for an infinite domain. In order to simulate finite-size problems, two approaches have been commonly used.
On the one hand, dislocation theory has usually tackled this by introducing image fields [18]. This approach can be difficult for general geometries, and it is rather complicated once the time variable has been introduced. On the other hand, the linear superposition scheme has been favoured by DDP. Its main advantage is the simplicity with which general geometries can be tackled, and the ability to solve the associated boundary-value problem using well-established numerical techniques such as the finite-element method. In this article, the superposition method shall be used as well.

(i) The linear superposition scheme
The application of the superposition principle to dislocation dynamics was first proposed by Lubarda et al. [34] and used widely thereafter by researchers following the Needleman & Van der Giessen [9] approach to DDP. The principle is summarized in figure 5.    The superposition principle is exact for each instant in time, so it requires no modification in the dynamic case described here. Let Ω be the domain of the boundary-value problem and Γ ≡ ∂Ω be the boundary of Ω, and let σ (x, t), u(x, t) be the stress and displacement fields therein. By virtue of the linear superposition principle, these can be conceived as the sum of two fields: σ =σ +σ and u =ũ +û, whereσ andũ are the stress and displacement fields of the infinite domainΩ containing the dislocations, andσ andû are the stress and displacement fields of a finite-size, dislocation-free mediumΩ coinciding with Ω. SurfaceΓ inΩ coincides with Γ and, as a result of the distribution of dislocations, there will be a tractionT and a displacementũ.
In order for the superposition of the fields inΩ andΩ to equate the fields in Ω, theΓ surface must have, for every instant in time, −T and −ũ applied over it in addition to the boundary conditions applied on Γ . Thus, if the dislocation-free domainΩ experiences both the actual boundary conditions and the reversed tractions and displacements caused by the infinite domain dislocation fields, the finite-size problem can be tackled as usual, with no image fields being necessary.
Furthermore, the elastic fields inΩ can be obtained through linear superposition of each dislocation's infinite domain fields, whereũ i ,σ i and˜ i denote, respectively, the displacement, stress and strain fields of the ith dislocation, N is the total number of dislocations andũ,σ and˜ are the total displacement, stress and strain.

(ii) The integration scheme
In general, the elastic fields in theΩ domain can be found using an elastodynamic numerical scheme such as finite-element or boundary-element methods; the elastic fields of theΩ domain are provided by the formulation presented here. However, one cannot describe the evolution of the dislocation structure over time without the addition of a mobility law and interaction rules. D3P can achieve a full description of the evolution and interaction of dislocation structures in two dimensions by adopting a forward Euler scheme as follows.
(i) At time t, for a given dislocation structure, calculate the total stress field σ =σ +σ = iσ i +σ . (ii) Calculate the Peach-Koehler force acting on all dislocations, where, for the ith dislocation, σ −σ i is used in order to omit dislocation self-stress; the same is done for all sources and obstacles using σ . (iii) Resolve the creation, annihilation, motion and interaction of dislocations with each other and obstacles according to short-range constitutive rules, updating their positions according to the mobility law for a time step t. (iv) CalculateT andũ onΓ for the updated dislocation structure.

Implementation rules and potential issues (a) Evaluation of the mobile contributions
The evaluation procedure for the elastic fields of the injection contributions does not differ from that of the quasi-static case. However, the mobile contributions are expressed as time derivatives of integral expressions that depend on the inverse of the past history of the dislocation, η(x), and therefore are of a different qualitative nature.  (x, 0). It is the presence of η(x) that makes the evaluation of the mobile contributions so fundamentally different from the static contributions: it implies that the elastic fields of the dislocations at some point in space and instant in time will depend not only on the current position of the dislocation, but also on each past position. Furthermore, because the elastic fields do not travel instantaneously, all the interactions of dislocations with one another and with the external fields will be based on a retardation principle: dislocations will not interact with the current position of others, but with their past history up to a certain past instant in time. Brock [38] provided a graphic representation of the influence of the past history function over the dislocation's surrounding medium. Pillon et al. [44], using a Peierls-Nabarro-Galerkin model, have addressed past history and retardation effects as well.

(i) The integration limits
The past history dependence and retardation effects are more easily appreciated by considering the evaluation of the mobile contributions. Consider the following integral: T ar 6 dξ , (4.1) Despite the integration limits suggesting otherwise, causality is observed: the integrand is multiplied by a Heaviside function, H(t − ar), that cancels the integrand for those values of ξ that, for a given spatial point (x, z) and instant in time t, the elastodynamic perturbations cannot have reached. Thus,t −ra plays the role that retarded times play in electrodynamics.
The integral above can therefore be thought of as z[t 2 (z 2 − 3x 2 ) + a 2 (2x 4 +x 2 z 2 − z 4 )] T ar 6 dξ , (4.2) where the integration limits ξ 0 and ξ t are functions of both t and (x, z) such that they cancel the retarded time. That is, ξ 0 , ξ t are such that It is clear that ξ 0 = 0 for all subsonic motion. As for ξ t , it marks the position the earliest radiated elastic perturbation might have reached at time t.

(ii) The past history function
For a generalized motion such as the one expected in a D3P simulation, the past history function η(ξ ) cannot be known a priori, so an analytic expression of the mobile contributions cannot be achieved in general.
Nevertheless, there are features of the past history function that are already known, at least as they would be obtained in any DDD simulation. In DDP simulations, dislocation motion is discretized into time steps, using the forward Euler approach [9] outlined in §3f. At a given time t 0 , the dislocation line will be at position x 0 ; after a time step of magnitude t, the current time will be t 1 = t 0 + t, and the position of the dislocation will have been updated to x 1 = x 0 + x 1 . Because the motion is discretized in time and space, interpolation between successive steps should respect whether the mobility law is linear or nonlinear. That is, η(ξ ) should have a form such that both the dislocation position and its temporal derivatives (velocity, acceleration, etc.) fulfil the mobility law over each time step.

(b) Numerical integration schemes
Because analytical expressions for the primitives of integrals such as that shown in equation (4.2) are not attainable except for relatively simple forms of η(ξ ) such as uniform motion [33], a numerical integration scheme is required.
It is most convenient and numerically effective for the integration interval (0, ξ t ) to be subdivided into the discrete steps of the past history function, and numerical integration performed within each of these steps. That is, where each (ξ i , ξ i+1 ) is one of the D3P integration intervals. The numerical method for the solution of each of the integrals is a matter of choice; in this work, an adaptive quadrature method based on Simpson's rule has been used. However, no matter which quadrature scheme is followed, special care should be taken due to the presence of several singularities in the integration path. Furthermore, there are two additional potential issues that must be addressed.

(c) Integration of the stress fields of the mobile contributions
As can be seen in table 2, the stress component fields of the mobile contributions are expressed as temporal derivatives of integral expressions. Direct numerical differentiation of these integrals is expensive if an acceptable degree of accuracy is required. Moreover, the solution consists of a series of radiated wave pulses, so the use of standard numerical differentiation schemes based on spatial (or temporal) discretization will be problematic. Furthermore, as Markenscoff & Clifton [33] argue, the interchange in the order of integration and differentiation is not, in general, legitimate. It would give rise to terms of the order of T −3 a and T −3 b that are, in principle, non-integrable singularities [33].
As Markenscoff & Clifton [33] proposed, consider a troublesome term taken from the transverse wave component of σ xz , There are two separate contributions. The second one does not produce a T −3 b singularity, so the order of integration and differentiation can be interchanged directly, In turn, as proposed by Markenscoff & Clifton, the first term can be integrated by parts [33], so that the resulting integral also allows the interchange with T b 0 = t 2 − b 2 r 2 , r = x 2 + z 2 and (4.8) where It must be pointed out that, if the motion of the dislocation is quiescent at t < 0, then η (0 + ) = ∞, so terms such as should vanish. A dislocation jumping from rest to some velocity would, on the other hand, have η (0 + ) = ∞, although this case has little physical significance. This procedure can be equally applied to the remaining terms, and it is necessary for a successful evaluation of the mobile contributions, irrespective of the numerical solution procedure chosen.
The special case of z = 0. The z = 0 case requires further attention. Consider, for instance, the following integral appearing in σ xz : (4.10) For z = 0, The integrand is singular ifx = 0, i.e. if x = ξ . This is not problematic as long as the singularity does not fall within the integration path; however, for points in x lying ahead of the injection site (x = 0) but behind the current position of the dislocation line, x − ξ is zero within the integration path.
Markenscoff & Clifton [33] and Markenscoff [37] found a way around this special case. Consider the example given by equation (4.10). This expression is derived from the following integral in Laplace space: The first inverse Laplace transform in space leads tô The usual way forward relies on interchanging the order of integration from ξ to λ; in that way, a single integral in its Cagniard form can be achieved. The interchange of integration variable can be performed only if Fubini's theorem is valid, but this theorem is applicable only if both the integrals in ξ and in λ converge [45]. Equation (4.10) highlights that the latter is not true for the aforementioned positions of x.
In order to perform the inversion, Markenscoff [37] proposed adding and subtracting the following term in Laplace space (for z = 0): so that the first reverse inversion leads to which amounts to adding and subtracting a term depending on the first-order Taylor expansion of η(ξ ) about ξ = x. By subtracting this term from equation (4.13), the singularity is cancelled; the term that is added corresponds to that of a uniformly moving dislocation, 2 and this problem was in fact solved analytically by Markenscoff & Clifton [33]. This solution leads to even more protracted formulae. In the example above, the solution would be and T a = t 2 − a 2x2 . Note that once this is done, I(x, 0, t) must be differentiated in time in order to obtain the transverse part of σ xz (x, 0, t). If necessary, this will be done by integrating by parts as already explained in §4c. Thus, Note that the z = 0 case is relevant for the σ xz and u z components alone; σ xx , σ zz , as well as u x vanish for z = 0.

(d) Singularities at the injection front and behind the injection front
The divergence of the elastic stress field at the dislocation core is a well-known feature in dislocation theory; as atomistic simulations show, it lacks any physical interpretation apart from the failure of first-order elasticity to determine faithfully stress fields close to geometric discontinuities due to the simplifying linearizations employed [46]. However, in computer simulations of dislocation dynamics, the presence of an infinity is not permissible, and several approximations have been proposed to deal with these singularities. In DDP simulations, the singularities at the core are simply removed [9]: a radial cut-off a few Burgers vectors wide is defined around the core, and the stress field within is assumed to be constant and equal to that on the boundary of the cut-off (taken radially).
As discussed in more detail in the electronic supplementary material appendix 'Locations of the singularities at the front and behind the front', in this case, there are two types of singularities: those due to the dislocation core and those due to the propagating fronts (cf. [47]). In order to prevent infinities from arising in DDP, a typical cut-off radius of 2-10 Burgers vectors is most commonly defined [9]. It is proposed that a D3P simulation should enforce a cut-off radius of similar magnitude around the moving core.
A safety ring also needs to be established around the singular regions of the injection front. If θ is the angular polar coordinate centred around the core, the locations over which the front's singularities vanish can be proved to be, for the longitudinal fronts: for σ zz , θ = 0, π ; for σ xx , θ = nπ; for σ xz , θ = (2n + 1)π/2. Here, n ∈ Z. Equally, for the transverse fronts, for σ zz , θ = nπ and θ = (2n + 1)π/4; for σ xx , θ = nπ and θ = (2n + 1)π/4; for σ xz , θ = (2n + 1)π/2 and θ = (2n + 1)π/4. This angular dependence can be appreciated in figure 6. Other than in these locations, where the fronts are non-singular, it is proposed that, as with the core cut-off, a ring of 2-10 Burgers vectors radius should suffice. Furthermore, the presence of such singularities at the front implies that if it encounters another dislocation, this dislocation will briefly be subjected to a very large Peach-Koehler force. Under such a large force, the mobility law must prevent the dislocation from reaching unphysical velocities.

Numerical results for the elastic fields of injected dislocations
The previous sections have outlined the derivation and method of evaluating the elastodynamic fields of an injected, non-uniformly moving straight edge dislocation. In this section, a number of illustrative numerical tests of this formulation are presented.

(a) The injected static dislocation
The stress fields shown in figure 6 display the characteristics that can be anticipated from their analytic expression: two independent monochromatic wavefronts are propagating outwards from the injection site at the longitudinal and transverse speeds of sound. In contrast to the elastostatic 'steady-state' solution (as derived in §3c), the elastodynamic fields do not exist everywhere in the domain. Once the transverse wave has passed a point, the magnitude of the field is seen to approach rapidly the values of the elastostatic solution at that point. This is seen clearly in figure 7, which shows the temporal evolution of the σ xz (x, z, t) component of stress: spreading outwards from the core, the magnitude of the field in figure 7c converges to that of the elastostatic solution in figure 7d after a few nanoseconds. However, the elastostatic solution becomes recognizable only after the transverse front has passed.

(b) The non-uniformly moving injected dislocation
The solution of the non-uniformly moving, injected edge dislocation consists of the superposition of the injected and mobile contributions. The former produces a set of two radiating fronts centred at the injection site and expanding outwards; as noted by Brock [38], the singularity at the transverse wavefront of the injection contribution cannot be cancelled by the mobile contribution  because the former is being superimposed on the longitudinal contributions from the radiating mobile core. Figure 8 shows the σ xz stress component of a uniformly and a non-uniformly moving injected dislocation. The uniformly moving dislocation advanced at 1000 m s −1 , whereas the speed of the non-uniformly moving one was allowed to vary randomly per time step in between M t = 0 and M t = 1. Here, M t = v/c t is the transverse Mach number, the ratio between the transverse elastic and inertia forces in the system. As can be observed, the impact of the past history over the values of the elastic fields is great, the randomized field displaying a great amount of irregularities with respect to the uniform one. Furthermore, the solutions in figure 9 show the σ xx and σ zz stress components for a random motion between M t = 0 and M t = 0.62 to ensure an average speed of v = 1000 m s −1 . Irregularities in the fields were reduced as the maximum speed was well below the transverse sound barrier.

(i) Peach-Koehler forces
Mura [4] proved that the Peach-Koehler force, which is usually derived in an elastostatic framework, remains valid when the dislocation is moving. The Peach-Koehler force exerted by one dislocation on another is directly proportional to the resolved shear stress of the former on the slip plane and in the direction of the Burgers vector of the latter. The time dependence of the force exerted by one dislocation on the other shows the same time dependence as the appropriate resolved shear stress of the former. Figure 10a shows the force at three instants in time exerted by an injected static dislocation at the origin on a second dislocation, the position of which is indicated on the horizontal axis. The dislocations share the same slip plane (i.e. z = 0), but they have Burgers vectors of opposite sign. The normal wavefront does not produce a singularity, but the transverse (shear) front certainly does. The intensity of the peaks of force ahead of the latter decreases with separation from the injection site.   Figure 10b shows the force exerted on a second dislocation, the position of which is indicated on the horizontal axis, by a dislocation that is moving with a constant velocity v glide = 1000 m s −1 . Again, the two dislocations share the same slip plane, but their Burgers vectors have opposite signs. The shear stress field of the moving dislocation is similar to that of the static solution, but the movement leads an extension of the field behind the dislocation compared with that ahead of it. Two solutions are displayed: one, based on the analytic solution given by Markenscoff & Clifton [33]; the second is the numerical solution achieved using an adaptive Simpson's rule quadrature of equation (3.39) for z = 0.

(c) Annihilation by an injected static dislocation
Consider the following situation: an injected moving dislocation is annihilated by the injection of a static dislocation of equal and opposite sign at its current position. What follows after some time is represented in figure 11: the injection of the static dislocation brings to an end the movement of the first dislocation, and waves are emitted from the position of the annihilation that cancel identically those of the moving one. In figure 11, one can observe four fronts: the outermost one  corresponds to the longitudinal front of the moving dislocation; the second outermost one is the longitudinal front of the opposite signed static dislocation; the third one, the transverse injection front of the moving dislocation; the innermost one, the transverse injection front of the static dislocation. The fields within that last circumference vanish. In fact, it can be analytically shown that in the limit t → ∞, their fields cancel, as expected. However, relevant to D3P simulations is that, as suggested with figure 11, once two unlike signed dislocations meet their annihilation, they will not immediately cancel the effect of each other's past history within the domain.

Conclusions
This paper introduces an elastodynamic formulation of DDP in which the time dependence of elastic fields of injected and moving dislocations is treated explicitly. This is a new paradigm in the simulation of DDD. The name dynamic discrete dislocation plasticity, or D3P for brevity, has been coined for this new formulation.
It has been shown that the application of standard quasi-static DDP to shock loading leads to the unphysical artefact of dislocation sources being activated ahead of the shock front. The reason for this failure is that quasi-static methods ignore the finite time it takes for elastic signals to travel in the medium: stresses created by dislocations behind the shock front are felt instantaneously by sources ahead of the shock front.
The two-dimensional dynamical formulation that has been described treats edge dislocations in plane strain. The solution for an existing non-uniformly moving edge dislocation, developed by Markenscoff & Clifton [33], has been augmented with the solution for an injected static edge dislocation. These two fundamental solutions are sufficient to model dynamically the generation, annihilation and movement of edge dislocations at high strain rates. The complexity of this formulation is considerable, especially in its numerical implementation, the subtleties of which have been discussed in detail; in particular, some of the difficulties that arise with singularities in the elastic fields and how they may be treated have been described.
The motivation for this work was to simulate the behaviour of dislocations in materials subjected to strain rates of 10 6 s −1 or more. Such high strain rates occur in shock loading. They may also occur in twinning and martensitic transformations, dynamical fracture, laser shock peening and in radiation damage by high-energy neutrons.
One of the key uncertainties in this work is the relationship between the damping force on a dislocation and its velocity, equation (1.2). This uncertainty arises in a variety of other contexts, such as the initiation and growth of fatigue cracks. It is a topic worthy of more research.