A One-Dimensional Peridynamic Model of Defect Propagation
and its Relation to Certain Other Continuum Models
The peridynamic model of a solid does not involve spatial gradients of the displacement field and is therefore well suited for studying defect propagation. Here, bond-based peridynamic theory is used to study the equilibrium and steady propagation of a lattice defect – a kink – in one dimension. The material transforms locally, from one state to another, as the kink passes through. The kink is in equilibrium if the applied force is less than a certain critical value that is calculated, and propagates if it exceeds that value. The kinetic relation giving the propagation speed as a function of the applied force is also derived.
In addition, it is shown that the dynamical solutions of certain differential-equation-based models of a continuum are the same as those of the peridynamic model provided the micromodulus function is chosen suitably. A formula for calculating the micromodulus function of the equivalent peridynamic model is derived and illustrated. This ability to replace a differential-equation-based model with a peridynamic one may prove useful when numerically studying more complicated problems such as those involving multiple and interacting defects.
Keywords: Peridynamic theory, defect propagation, kink, phase transformation, dynamics, Frenkel-Kontorova.
In a seminal paper , Silling introduced a nonlocal theory of elasticity in integral form that accounts for the interaction between continuum particles at finite distances, reminiscent of models involving interatomic interactions. The theory, called Peridynamics, has subsequently been extended to situations involving non-central forces and inelasticity; see for example the review article . It has recently gained significant attention, at least in part because of the advantages it affords numerical calculations involving discontinuities111For example, Sun and Sundararaghavan  show in their study of crystal plasticity that peridynamics can model finer shear bands with less simulation time than the traditional finite element solution.. Since the theory does not involve displacement gradients, it can handle geometric singularities with relative ease and is therefore particularly well suited for the study of defects. This is the motivation for the present study where we analytically examine the propagation of a certain type of defect – a kink – according to peridynamics.
Analytical solutions within the peridynamic theory are few, and limited to linear problems, such as determining Green’s functions, Weckner et al. , Wang et al. , and using them to study the static and dynamic response of an infinite bar, Silling et al. , Weckner and Abeyaratne  and Mikata . While the peridynamic operator in our study is linear, the effective body force is nonlinear, because accounting for the defect requires the material to have two states, and therefore the model to involve a double-well potential. However, by taking this potential to be bi-quadratic, we are able to analytically solve the peridynamic defect propagation problem.
Energy functions with multiple local minima (energy-wells) are frequently encountered when modeling various physical systems. In equilibrium, the system maybe “stuck” in a metastable energy-well, and under a suitably large external stimulus, e.g. force, will undergo a progressive transition from the metastable state towards a stable state. Often, such an evolution involves a kink transition front that transforms the system locally as it passes through. Such phenomena are observed in, for example, mechanical systems (friction , dynamics of CNT foams , lattices of bistable buckled elastic structures , and mechanical transmission lines ); material systems (dislocation motion , ferromagnetic domain wall motion , commensurate phase transitions , and chemical surface adsorption ); electromagnetic systems (magnetic flux propagation in Josephson junctions ); biological systems (pulse propagation in neurophysiology  and rotation of DNA bases ); and even traffic flow . Determining the speed with which the system transitions into the stable state, as a function of the applied stimulus, is a question of significant interest.
There is a considerable literature on the propagation of kink transition fronts in lattices. The model proposed by Frenkel and Kontorova (FK) , that now goes by their names, underlies many such studies. Atkinson and Cabrera  made an important contribution in the analysis of the FK model where, by using a bi-quadratic double well potential (rather than a trigonometric one) they solved the problem in analytic closed form. In particular, they derived a relation between applied force and kink propagation speed, and demonstrated the apparent dissipation of energy due to its radiation by waves propagating away from the kink. Several important contributions concerning kink propagation through a discrete lattice have been made by Truskinovsky, Vainchtein and their collaborators, e.g. [19, 20, 35, 36, 37, 38, 39]; for a recent study, see Shiroky and Gendelman . A survey of this literature can be found in the book by Braun and Kivshar ; see also Slepyan . A striking characteristic of Atkinson and Cabrera’s solution is that it predicts unbounded values of force at certain small values of speed.
Turning next to continuum models of kink propagation, the simplest is where the lattice is replaced by a linear elastic solid (but with a bi-quadratic double-well onsite potential). Here one finds that the force must equal the Maxwell force for all subsonic propagation speeds, and must take a different value for all supersonic speeds, e.g. see [2, 13, 14]. This indeterminacy of the speed corresponding to a given value of force is likely due to the absence of a length scale in classical elasticity; after all, the discreteness of a lattice model introduces a length scale into the problem even if only nearest neighbor interactions are accounted for. The next simplest continuum model is obtained by replacing the lattice by an elastic solid that includes strain gradient effects – the so-called Boussinesq approximation. Kink propagation according to this model has been studied in [2, 13, 14] but, as noted by Kresse and Truskinovsy , is a poor approximation at small wave lengths due to instability. Kresse and Truskinovsy  and Truskinovsky and Vainchtein  have studied an alternative model that does not suffer from this deficiency, the so-called quasi-continuum model that is based on an approach put forward by Rosenau . The quasi-continuum model may be viewed as arising from the retention of ‘micro-inertia” terms in Mindlin’s strain gradient theory . More elaborate kink propagation models based on more detailed atomistic considerations have also been studied, e.g. Abeyaratne and Vedantam  and Hildebrand and Abeyaratne . Results pertaining to FK models in higher dimensions are described in Chapter 11 of .
In the present paper a discrete lattice is approximated by a peridynamic model, and this is used to study kink propagation. The peridynamic operator is taken to be linear, bond-based and involves a single “micro-modulus function”. However the governing equation is nonlinear due to the presence of the double-well onsite potential that occurs in the body force term outside the peridynamic operator222Double-well potentials arise in modeling both kink and phase boundary propagation, the mathematical distinction being that in the latter case the nonlinear potential occurs within the peridynamic operator as in Dayal and Bhattacharya . The analysis in this case is considerably more difficult, e.g. the solution in  is numerical in contrast to the analytical solution we find here.. Motivated by Atkinson and Cabrera , the double well potential is taken to be bi-quadratic. This introduces three parameters: a modulus , distance between energy wells , and the height of the energy barrier . As for the micromodulus function, its simplest form involves two material parameters, an elastic modulus and a length scale , whose values are chosen by matching features of the dispersion relations of the peridynamic and discrete lattice models. The particular form of the micromodulus function we use is motivated by Silling , and remarkably, is also implied by the quasi-continuum model as we shall show. Using Fourier transformation, we solve in closed form the static problem for an equilibrium kink and the dynamic problem for a steadily propagating kink. A kink will be in equilibrium if the applied force is smaller than a certain value333The fact that exceeds the Maxwell force is said to describe the phenomenon of “lattice trapping”. . When the force exceeds this value, the kink can propagate at a steady speed given by a kinetic law that is determined by the analysis. This law relates the applied force to the propagation speed. As the speed increases, the force needed to propagate the kink increases monotonically from to . The dynamic solution involves effective damping caused by energy radiation away from the kink and is well-behaved at all velocities.
After solving the kink propagation problem using peridynamic theory, we observed that the results had a similar form to those arising from the quasi-continuum model. Indeed, we found that for a particular choice of the ratio of the peridynamic and lattice length scales, , the results coincided exactly. This led us in Section 5 to look at the relation between peridynamic theory and differential equation models of elastic continua more generally in settings not limited to kink propagation. In Section 5 we do not assume the form of the micromodulus function a priori, but rather determine it such that the results of the peridynamic theory and a class of differential-equation-based theory coincide. Using an approach based on Whitham , a formula for calculating the equivalent micromodulus function is derived. It is then possible to replace a differential-equation-based model with an equivalent peridynamic one, and this may prove useful when numerically studying more complicated problems such as those involving multiple and interacting defects.
The paper is organized as follows: first, the basic one-dimensional discrete model of kink propagation through a lattice is described in Section 2. For purposes of comparison, in Section 3 we state (without proof) the results of the classical linear elastic solution of this problem. In Section 4 we describe the linear bond-based peridynamic model; show how its parameters can be chosen by comparison with the lattice model; and then go on the solve the static and dynamic problems for a kink. In particular, the kinetic law relating the applied force to the kink speed is derived. In Section 5 we compare the results of the previous section with those of the quasi-continuum theory, and this leads us to explore the relation between the peridynamic model and other differential-equation-based continuum models in a more general dynamical setting.
Consider a one-dimensional row of particles. In a reference configuration, they are equally spaced, a distance apart, with the th particle located at , . During a motion, the displacement of the th particle at time is .
Each particle has mass and is connected to its nearest neighbors by linear springs. Specifically, for each , the th particle is connected to the th particle by a spring of stiffness and unstretched length . Each particle is also attached to a fixed foundation by a nonlinear spring characterized by a double-well potential as shown schematically in Figure 1. The local minima of are at and . The effect of the external loading is accounted for by a force acting on each particle. The equation of motion of the particle is therefore
It is convenient to write
where is a suitably smooth function defined for . Then (1) can be written as the delay-differential equation
where we have let and .
In order to obtain explicit closed-form solutions to the problem to be studied, we take the double-well potential to be bi-quadratic:
where and are constant parameters. Observe that has local minima – energy-wells – at and , and a local maximum at . If the displacement of some particle takes a value we say it is in “phase-1”; if it exceeds , we say it is in “phase-2”. The force444For simplicity of terminology we refer to and as “forces” though they are in fact force densities: forces per unit volume. associated with this potential is bi-linear:
It is convenient to let and denote the force levels
whose meanings are shown in Figure 2. A particle cannot be in phase-1 if the force exceeds , and it cannot be in phase-2 if the force is smaller than . The Maxwell force-level cuts off lobes of equal area as shown in the figure. It is the force at which the two phases have the same potential energy.
In the particular case of a uniform equilibrium configuration in which all particles have the same displacement , it follows from (2) that is given by the roots of the equation . For the bi-quadratic energy (3), when the applied force lies in the intermediate range
(see Figure 2), this leads to two possible equilibrium displacements :
In Sections 3 and 4 our primary interest will be in the steady motion of a “kink” associated with the motion of particles from phase-1 to phase-2, one at a time. At a typical instant, all particles ahead of the th one are still in phase-1 while those behind it have moved into phase-2. Suppose that during this motion the particle stays in phase-1 for times (for some ); at time it moves into phase-2 and remains there for :
This motion continues steadily with the th particle moving into phase-2 at time and so on. Such a steady motion is described by
Setting in (8) shows that where is the propagation speed of the kink. Therefore the displacement field of interest has the traveling wave form
The particles far ahead of the kink are taken to be in equibrium is phase-1 so that as , i.e. as .
We use the term “kink” to refer to a configuration, either static or dynamic, in which all particles ahead of a particular one are in phase-1 and those behind it are in phase-2 .
3 Classical elasticity model.
The simpest mathematical model of the problem described above is obtained by replacing the finite difference terms in (2) by the second spatial derivative of leading to
The associated material wave speed is Steady traveling wave solutions, , of this equation have been studied, e.g. see [2, 13], with particles ahead of being in phase-1, those behind it in phase-2, and
The displacement and strain, and , are required to be continuous. In the case of subsonic propagation the solution is
where (10) follows from the requirement necessitated by the continuity of . Thus in the subsonic case, no matter what the propagation speed, the force must equal the Maxwell force .
On the other hand for supersonic propagation speeds one finds
Thus in the supersonic case, no matter what the propagation speed, the force must equal the maximum force level indicated in Figure 2.
The horizontal solid lines in Figure 3 show the relation between the force and propagation speed as described by (10) and (12). (The curve corresponds to the peridynamic solution to be determined later.) Observe how in the subsonic case the solution (9) decays exponentially away from the kink in both directions, but in the supersonic case involves oscillatory waves that are radiated without decay behind the kink, see (11). The radiated energy in these waves describes “dissipation”; this is also implied by the fact that exceeds the Maxwell-force in this case.
4 Peridynamic model.
A peridynamic model is a nonlocal continuum model that accounts for the effect of a neighborhood on a material point. In (linear, homogeneous, bond-based) peridynamics, the force per unit length applied on the particle at by the particle at is
and so the equation of motion (2) is replaced by its peridynamic counterpart
The micromodulus function in (13) characterizes the material. On physical grounds is required to be symmetric: ; to decay at infinity: as ; and to involve a length scale such that approaches a suitable generalized function in the limit .
Note that since particles interact at a distance, it is not necessary that the displacement field in peridynamic theory be continuous. For example suppose that the micromodulus function is nonzero for and vanishes identically for for some (the “horizon”). The continuum would then remain coherent even in the presence of displacement discontinuities of magnitude less than .
Throughout this section we shall consider a micromodulus function of the form
where the longitudinal stiffness and the length are material constants. This form is motivated by an analysis of nonlocal effects by Silling . It will be convenient in what follows to let Deleted from next equation.
and the intrinsic wave speed associated with this material model by
Observe that if instead of (14), the micromodulus function is taken to be
where denotes the Dirac -function, the peridynamic equation of motion (13) coincides with the equation of motion (2) of the discrete model. If we want (14) to approximate (17), we can choose the material parameters and in (14) by matching, for example, the dispersion relations of the two models. By looking at motions of the form of equations (2) and (13) (with and set equal to zero) one obtains the following respective dispersion relations relating the frequency to the wave number :
In the shortwave length limit, , the peridynamic dispersion relation (18) gives
The discrete dispersion relation oscillates in this limit and so the value of can be chosen by matching some desirable average feature. Figure 4 shows plots of the dispersion relations (18) and (18). In drawing the figures we have taken , , and plotted versus . Three peridynamic dispersion relations are shown corresponding to three values of illustrating how one can vary the way in the peridynamic and discrete dispersion relations match by varying this parameter. One could of course work with a micromodulus function with additional parameters and match additional details of the dispersion relations.
Finally, it is useful to introduce the function , defined in terms of the micromodulus function by
where the second and third equalities follow since is an even function. will play a central role in what follows (and already appeared in (18) above). It is a one-dimensional scalar counterpart of the acoustic tensor, see , and so we will refer to it as the acoustic function. When the micromodulus function has the exponential form (14), the acoustic function specializes to
4.1 Equilibrium configuration of a kink.
For an equilibrium configuration , the peridynamic equation of motion (13) specializes to the equilibrium equation
We now consider non-uniform equilibrium configurations in which all particles are associated with phase-1, and all particles are associated with phase-2, with
For the class of displacement fields under consideration,
Observe from (20) that because of the discontinuity of , must be discontinuous at and so we are not requiring in (22). Writing (22) as where is the Heaviside step function, and using this in the bi-quadratic energy function (3) leads to
Therefore the equilibrium equation (20) specializes to
As observed already, the discontinuity of at requires to be discontinuous at . To calculate the jump in the value of at , we evaluate (23) at and at and subtract one equation from the other. Letting in the result leads to
where we have used the standard notation for the jump in the value of a function at :
Thus the displacement suffers a jump discontinuity at of magnitude
where was introduced in (15). Displacement jumps in the peridynamic theory were encountered in , where it was shown in particular, that such discontinuities cannot propagate. Indeed, when we consider the propagating kink in the next section, we find no displacement discontinuities.
and we have set
The integral in (27), interpreted as its Cauchy principal value, can be evaluated by contour integration using the Residue Theorem.
To do this, note first by (28), that the integrand in (27) has three poles, one corresponding to and the other two being the two (purely imaginary) zeros of . Next, one considers a path along the real axis of the complex -plane, indented by a small semi-circle in the upper half plane555Indenting in the upper half-plane rather than the lower half plane is necessary in order to get the correct conditions at infinity. of radius centered at ; the complex -plane and integration path are shown in Figure 6 of the Supplementary Material. The desired integral is then the limit as of the difference between the integral on minus the integral on . The integral on , in the limit , is readily shown to be . Thus (27) can be written as
This integral on can now be evaluated by completing the integration path by a large semi-circle in the lower half-plane for and the upper half plane for and using the Residue Theorem. This leads to
where is the Maxwell force given in (5). It is not difficult to show using (5) that the rightmost expression in (32) is while the leftmost expression is . Thus the interval of force demarcated by (32) is a subset of the interval that includes the Maxwell force in its interior.
4.2 Dynamics of a steadily propagating kink.
We now consider the steady motion of a kink propagating at some (to-be-determined) speed . The displacement field in such a motion has the form
Substituting this into the peridynamic equation of motion (13) gives
The material ahead of the propagating kink is in phase-1, that behind it is in phase-2:
According to (33), in the dynamic case the discontinuity in at requires that be discontinuous at but and are permitted to be continuous. Therefore, in particular, we have required in (34). This is in contrast to the equilibrium problem. As for the jump in , it is seen from (33) that whence
Finally, as for the far-field conditions we assume that the particles far ahead of the kink are in equilibrium in phase-1. Thus, recalling (7), we take
As for the particles far behind the wavefront, we simply require them to be in phase-2 but do not assume them to be in equilibrium. Thus as we do not require anything beyond (34). This allows for particles to be oscillating in the phase-2 energy-well. We could impose the stronger requirement that the average displacement of the particles behind the wave be (where is given by (7)) but it is not necessary that we do so. The results show that this comes out automatically.
with and being the (positive real) quantities
Here we have let
where the integral is interpreted in the sense of its Cauchy principal value.
The integral in (42) can again be evaluated by contour integration using the Residue Theorem in a manner similar to that used in the study of the equilibrium kink in Section 4.1. The main difference is that here, the integrand has two poles on the real -axis (in addition to the pair of purely imaginary poles). When constructing the integration path we must indent it so that passes below these two poles. The complex -plane and integration path are shown in Figure 7 of the Supplementary Material. The reason it must pass below (rather than above) these poles is because these poles lead to oscillatory terms and in the solution, which cannot exist as since those particles are assumed to be in equilibrium. On the other hand we are allowing the particles to vibrate as and so such oscillatory terms are admissible for . By closing the path of integration as in Section 4.1 and using the Residue Theorem we find
It can be readily verified that and are continuous at and that has the appropriate discontinuity (35). The requirement is yet to be enforced.
An alternative illuminating form of the solution is
Observe using (5) and (7) that the coefficient of the exponential terms vanishes if while the coefficient of the cosine term vanishes if . The displacement field written in the form (44) satisfies the requirements and . The requirement on in (35) remains to be enforced.
4.3 Kinetic relation.
While given by (43) is continuous at we are yet to enforce the requirement necessitated by (34). When this is enforced (43) leads to the following relation between the applied force and the kink propagation speed :
Given the value of the force , the kinetic relation (46) gives the value of the propagation speed666Alternatively the kinetic relation follows by enforcing the requirement (35) on the solution (44). . Observe that and therefore that .
When the lattice length scale is much smaller than the peridynamic length scale , we let in (46) to see that at each fixed . On the other hand when the peridynamic length scale is much smaller than the lattice length scale , letting shows that for and for ; see Figure 3. Note that when at fixed .
Finally, recall that a kink can be in equilibrium if the applied force lies in the range given in (32). The rightmost expression there, i.e. the maximum value of that allows for an equilibrium kink, is precisely . Thus if the value of the force is monotonically increased from zero, the kink will remain stationary for and start to propagate once the force exceeds . Since it follows that the kink will remain in equilibrium even if the applied force exceeds the Maxwell force provided it is less than – the phenomenon sometimes referred to as “lattice trapping”.
4.4 Dynamic solution when .
Recall that the displacement field in the dynamic solution (43), (46) is continuous at , but is discontinuous at in the static solution (30). In order to understand this difference, we now examine the dynamic solution at small values of the speed777One can see from (48) and (49) that, at any fixed negative value of , the limit of as does not exist in the usual sense. It must be interpreted in the sense of Young measures. . It is straightforward to show from (40) and (43) that for small the dynamic solution takes the form
to leading order, where from (40)
The small- dynamic solution (48), (49) would be identical to the static solution (30), (29) if not for the cosine term in (48). At fixed (no matter how small) given by (48) is continuous at . However as observe that the frequency of oscillation of the cosine term goes to infinity whereas its amplitude remains strictly positive and finite. Observe also that the coefficient of the cosine term term is precisely equal to the jump in displacement of the static solution as given in (26). Of course the static solution (30) holds for all values of the force in the range (32) whereas the dynamic solution for small holds only for .
5 Equivalency of peridynamic and certain other continuum models.
In the preceding section we analyzed the kink propagation problem within the peridynamic theory, a problem that has been studied previously within the quasi-continuum theory by Kresse and Truskinovsky  and Truskinovsky and Vainchtein . It turns out that the displacement field and kinetic relation according to that theory have the same forms as the respective equations (43) and (46) of the peridynamic theory. In fact, one can show that they coincide exactly when the peridynamic length scale has the particular value , indicating that the quasi-continuum theory can be viewed as a special case of the peridynamic theory for a particular value of the length scale, at least in the context of kink propagation. In general, the presence of the parameter in the peridynamic theory gives it an additional degree of freedom, and its value can be chosen, for example, to match features of the discrete dispersion relation as in Figure 4. Even so, the aforementioned coincidence of results when is rather striking.
We did not anticipate this result, and it prompted us to ask the following more general question: given a continuum model described by some differential equation, does there exist an “equivalent” peridynamic model? This question is stated more precisely and investigated in the present section. The discussion pertains to any dynamical motion on an infinite interval (not only steady kink propagation). Moreover, we will not pick a particular peridynamic micromodulus function a priori. Instead, we shall determine the function that makes the peridynamic model equivalent to the other model described by a differential equation. This will be illustrated using two explicit examples (the Boussinesq and quasi-continuum models) after a more general analysis.
We shall use the notation for the Fourier transform of a function . Since different authors use slightly different definitions of this, we note that the one we use is
Our strategy is based on Whitham : we first Fourier transform the peridynamic equation of motion. The result involves the acoustic function which can be determined from (19) if we know the micromodulus function . We invert this relation and obtain a formula for calculating when is known. Then, given some other continuum model described by a differential equation, if the Fourier transform of that equation of motion has the same form as the Fourier transformed peridynamic equation of motion for some we say the models are equivalent, and we can calculate the associated micromodulus function provided that a certain integral exists.
We start with the peridynamic equation of motion
where the body force may equal as in the preceding section. While one could consider more general body forces, e.g. , one of the attractive features of peridynamic theory is that it involves no spatial derivatives of and so including a -dependency in the body force would be counter to that spirit. Therefore we limit attention to body forces of the form . Fourier transforming (51) leads to
where the acoustic function is
Given the micromodulus function , (53) is an equation for calculating the acoustic function . Our immediate goal is to invert this relation to calculate in terms of .
Observe that for continuous displacement fields ,
Therefore, adding a Dirac -function to the micromodulus function does not change the peridynamic equation of motion (51). Likewise it does not change the acoustic function given by (53). Thus we will only determine to within a -function.
First, Dayal  has shown that if has a finite limit as then
Thus we can write (53) as
which says that is the Fourier transform of . Thus by using the inverse Fourier transform (50) we obtain
which gives the micromodulus function corresponding to the acoustic function .
The preceding result is not valid if as . Suppose that as for some integer , with
being a finite number. In general, the micromodulus function will then involve a generalized function and its derivatives as in Whitham . Even so a formal characterization of is possible as follows: define the function by
One can show by taking the Fourier transform of this equation and integrating by parts that
Given a continuum model based on some differential equation, if the Fourier transform of the associated equation of motion has the form (52) for some , we say the models are equivalent, and one can determine the micromodulus function of the equivalent peridynamic model through either (56) or (58), (59), provided the relevant integrals exist.
In the next subsection two explicit differential-equation-based continuum models stemming from the basic lattice equation of motion (2) will be examined. They will both be of the generic form
where the differential operators and are
for some constant ’s and ’s, and integers . Recalling that for any function
where the polynomials and are