Fourier-based numerical approximation of the Weertman equation for moving dislocations

Fourier-based numerical approximation of the Weertman equation for moving dislocations

Marc Josien École des Ponts and INRIA, 6 et 8 avenue Blaise Pascal, 77455 Marne-La-Vallée Cedex 2, France.    Yves-Patrick Pellegrini CEA, DAM, DIF, F-91297 Arpajon, France.    Frédéric Legoll École des Ponts and INRIA, 6 et 8 avenue Blaise Pascal, 77455 Marne-La-Vallée Cedex 2, France.    Claude Le Bris École des Ponts and INRIA, 6 et 8 avenue Blaise Pascal, 77455 Marne-La-Vallée Cedex 2, France.
July 5, 2019

This work discusses the numerical approximation of a nonlinear reaction-advection-diffusion equation, which is a dimensionless form of the Weertman equation. This equation models steadily-moving dislocations in materials science. It reduces to the celebrated Peierls-Nabarro equation when its advection term is set to zero. The approach rests on considering a time-dependent formulation, which admits the equation under study as its long-time limit. Introducing a Preconditioned Collocation Scheme based on Fourier transforms, the iterative numerical method presented solves the time-dependent problem, delivering at convergence the desired numerical solution to the Weertman equation. Although it rests on an explicit time-evolution scheme, the method allows for large time steps, and captures the solution in a robust manner. Numerical results illustrate the efficiency of the approach for several types of nonlinearities.

Weertman equation, Peierls-Nabarro equation, dislocations, Cauchy-type nonlinear integrodifferential equation, reaction-advection-diffusion equation, fractional Laplacian, preconditioned scheme.

I Introduction

This article addresses the numerical approximation of the following nonlinear integrodifferential equation, with Cauchy-type singular kernel:


where both the real-valued function and the scalar constant are the unknowns. The potential is a nonlinear bistable function of with (at least) two local minima at values and . The meaning of the subscript is explained below. The operator is linear, and defined in terms of the Hilbert transform KING09 () as


where denotes the principal value KANW04 () at . In the context of singular integral equations, the above kernel of the Hilbert transform is known as the Cauchy kernel. The operator , also denoted by some authors, is diffusive (Brezis, , p. 181). Another useful representation of (2) is (by integration by parts)


Let the Fourier transform in the continuum (FT) of a function be defined at wavemode as


One has (Gradshteyn, , p. 1118), whence . Thus, the non-local operator is symmetric and positive.

Equation (1) is a dimensionless form of the Weertman equation WEER69a (); WEER69b (); Rosakis () (simply referred to as ‘the Weertman equation’ in the following), which models straight dislocations traveling with steady velocity, thus generalizing the Peierls-Nabarro (PN) equation for static dislocations Nabarro46 ():


Dislocations are linear defects in crystals, the motion of which is responsible for the plasticity of metals Hirth (). Dislocation lines have a non-vanishing sectional area, i.e., they possess a ‘core’ of finite width. The derivative (the so-called dislocation density) of the unknown function in (1) describes the shape function of a flat finite-width dislocation on its glide plane, along the -direction. The core is the region of space where develops peaks. From a physical standpoint, the function represents a local relative material displacement discontinuity between the upper and lower half-spaces surrounding the glide plane on which moves the dislocation line; see, e.g., Hirth () for details. From a broader perspective, the function can be understood as a moving phase-transformation front between the states and (Fig. 1).

Figure 1: Typical shape of in Equation (1) when is a sinusoidal function.

In (1) the term accounts for the long-range elastic self-interactions that tend to spread the core. This repulsive interaction is counterbalanced by the nonlinear pull-back force , which binds together the upper and lower half-spaces, thus giving the dislocation core its finite width. Throughout this article, we consider that includes a constant externally applied loading (that is, where is an energy potential intrinsic to the material, tilted by the adjunction of a linear term ). Moreover, the moving dislocation is subjected to various drag mechanisms encoded into the term . As recalled below, the Weertman equation admits an analytical solution Rosakis () if is assumed sinusoidal. For more realistic potentials, a numerical approach is required.

Since its inception, the original PN model has been enriched in various directions. For instance, it most often requires being generalized to vector-valued to be quantitatively predictive SydowHartfordWahnstrom (); LuBulatovKioussis (); LuBulatovKioussis2 (); DENO04 (); DENO07 (). Also, the model has been extended to two dimensions of space, to study planar dislocation loops DENO04 (); DENO07 (); Xiang (); zhu2015 (). Quite generally, methods to compute the shape of static or moving cores encompass variational approaches and involve finite-element and/or phase-field-type implementations DENO04 (); DENO07 (); ZHAN15 (); Mianroodi (). Yet, in spite of such a wealth of enrichments of the PN model and its associated numerical methods of solution, the one-dimensional Weertman equation —a comparatively simpler extension— has not been investigated as thoroughly, while the specific problem of determining the allowed velocities of steadily-moving dislocations for general force laws remains an open question of major practical interest Rosakis (). For this reason, the present work focuses on solving the simplest, scalar, and one-dimensional case. A generalization to the vector case of the Weertman equation is the subject of ongoing work, and will be presented elsewhere YPP2 ().

It must be emphasized that, in the dimensionless form (1) of the equation, is not the physical velocity of the dislocation. The latter is deduced from in a post-processing step to be explained in YPP2 () (it is a direct application of Equations (46), (47) or (48) of YPP ()), which in the present scalar case is independent from the numerical task of solving the equation. Therefore the physical velocity is not further considered hereafter. Moreover, we stress that (1) applies only to subsonic motion, for which the coefficient of does not vanish in the original Weertman equation.

Numerical methods for solving integrodifferential equations such as (1) have been proposed by many authors. The method employed in Lejvcek () uses properties of the Hilbert transform to recast Equation (5) into a form amenable to fixed-point approaches. In Kurtz () the authors consider a simpler version of (1) on a bounded domain, in which the nonlinear term is replaced by some given -independent function , and where is also given. The solution is then obtained by means of a collocation method with quadratic interpolation. Those works make use of the expression of the operator in the direct space. More recently, Karlin et al. presented Movchan () a general iterative method for solving (5), based on the expression of in the Fourier space. Our work borrows from the latter reference. The interested reader can also refer to BUEN14 () for Fourier-based numerical schemes applied to the fractional Laplacian operator with .

The present article proposes a numerical method to approximate solutions to (1) in the case where is bistable. As in Movchan (), we build a dynamical system that admits (1) as its long-time limit, namely,


where , and is a regular initial data taking values in the interval such that


where are zeros of .

The iterative numerical approach introduced below uses (6) to approximate the solution of (1). It is immediate that if solves (1) and if we impose , then satisfies (6) for the initial data . It is proved in MathWeert () that if is bistable, then for any initial data with values in that satisfies (7) and for any continuous function , the solution of (6) converges to the solution of (1) in a sense explained in Sec. II.5 below. This leads to the following procedure:

  1. consider an initial condition such that (7) holds;

  2. approximate the solution to (6);

  3. while evaluating , choose such that the core of remains within the computational box and that converges to ;

  4. for sufficiently large so that has numerically converged, return as a numerical approximation to , and deduce the velocity .

As will be shown, this strategy proves efficient in cases of physical interest. However there might exist alternative strategies for solving (1).

Numerical approximation of (1) paves the way to investigating numerically the Dynamic PN equation YPP (), which generalizes the Weertman equation to transient regimes. Indeed, the initial conditions and long-time steady-state regimes of the Dynamic PN equation are solutions to the latter equation YPP1 (). However, we emphasize that (6) is only an algorithmic tool that has no relationship whatsoever with the actual dynamics of dislocations.

The article is organized as follows. In Sec. II, we briefly study the Weertman equation, and discuss both the uniqueness of its solution and its interpretation as the long-time limit of the dynamical system (6). We formally derive asymptotes of solutions to (1) and state identities about the velocity in general cases. Also, we explain how to choose in (6) to solve this equation in a comoving frame, namely, one which follows the translational motion of the core. An analytical solution to (1) that exists in a particular case is recalled. In Sec. III, we introduce our numerical representation for and discuss corresponding implementations of the diffusive operator and the advective operator , as well as methods to evaluate . Also, we make use of the asymptotic behavior when of the solution to (1) to circumvent the issue of the infinite domain of integration in (2). Once these fundamental elements have been introduced, we build in Sec. IV a Preconditioned Collocation Scheme (PCS) that applies to our problem, and justify this denomination. In Sec. V, we use this numerical approach on two test cases: one with a simple potential , for which the exact solution is known, and one with a more physically relevant potential . We also illustrate the robustness of our approach, concluding that the algorithm presented is unconditionally stable with respect to the time step . We empirically derive error scalings with respect to the parameters involved in the discretization. A concluding discussion closes the article, underlining some limitations of our approach, and proposing a few possible extensions. An Appendix is devoted to examining further one such extension.

Ii Some properties of the Weertman equation

This section is devoted to an overview of some important properties of the Weertman equation (1) and of the companion dynamical equation (6).

ii.1 Invariances

As the PN equation, Equation (1) is obviously invariant by translation. When invoking uniqueness of the solutions, we shall henceforth implicitly refer to ‘uniqueness up to arbitrary translations’. This invariance has consequences on the numerical solution, which usually undergoes an undesirable drift during calculations if no corrective action is undertaken. A special procedure is developed in Section II.6 below to eliminate this difficulty.

Moreover, Eq. (1) is invariant by reflection in the sense that if is a solution for boundary conditions , then is another solution for boundary conditions . Therefore, without loss of generality, we always assume throughout this article that .

ii.2 Existence and uniqueness of solutions to the Weertman equation.

There exists a unique solution to (1) when is a bistable nonlinearity. More precisely it can be shown rigorously (the proof relies on a recent result Gui ()) that for sufficiently regular, if are such that: (i) and ; (ii) any local minimizer of between and satisfies and , then there exists a unique velocity and a decreasing function satisfying (1), which is unique up to translation. Condition (i) means that is a ‘bistable potential’. A typical example of such , to be used in Sec. V.5, is drawn in Fig 2. It corresponds to

Figure 2: Camel-hump potential defined by (8), with parameters and .

Possible non-decreasing solutions to (1) that might exist as in the PN equation (Nabarro46, , Equation (16)), e.g., for , are disregarded in the present work. Hereafter, we always assume that obey the above conditions, and we term such values of at infinity consistent boundary conditions (CBCs).

ii.3 Asymptotic behavior and characteristic lengths

By letting , we formally deduce the leading-order asymptotic expansions


They are proved rigorously in MathWeert (). The key ingredient of the proof is the following asymptotic behavior of integrals with Cauchy kernel (see (Musk, , p. 267)):


It can be formally retrieved from the leading term of a ‘series expansion’ of (2) at large (note however that the integral involved in the next-to-leading term of such a formal expansion may not exist). Substituting expression (10) into (1), using the Taylor expansions at boundary values


and noting that vanishes like Gui (), we obtain (9).

Introducing the characteristic lengths


Equation (9) can be rewritten as


It will be argued in Sec. II.7 that represent typical scales of variation of the dislocation density on both sides of the solution. For further purposes, it is useful to introduce a mean characteristic length as


Finally, depending on the potential at hand, the solution can be either symmetric, in the sense that is an even function (an example will be given in Sec. II.7), or non-symmetric if the latter property does not hold.

ii.4 Velocity determination

There are many ways to determine the velocity associated with the solution to Equation (1). One possibility is to multiply (1) by a function chosen such that the integrals involved make sense, and to integrate over . This yields


In particular, the following choices of eliminate . With one gets


where the improper integral is understood as a principal value at infinity (Gradshteyn, , p. 252). Taking instead Gui (), one arrives at an expression easier to implement than (16), namely,


ii.5 Convergence towards solutions to the Weertman equation

The first author (M.J.) proves in MathWeert () that, under our working hypotheses, all the solutions to (6) converge towards the unique solution of (1) at exponential rate. Indeed, consider the following equation:


where is the initial condition of (6). The connection between (18) and (6) resides in that solves (18) if and only if


solves (6). Now, under mild requirements similar to those of Sec. II.2, Equation (18) can be shown MathWeert () to have a unique solution with the following property: if is the solution to (1) with same boundary conditions at infinity as , then there exist constants , and such that


The proof relies on a comparison principle, which is a generic property of operators where is a dissipative operator (see the classical references Fife (); XinfuChen ()).

Combining Equations (19) and (20), one deduces that at large times and uniformly in ,


Thus, given CBCs at infinity, the long-time limit of the solution to (6) is the solution to (1) with same CBCs, up to a time-dependent drift. In the next section, we show how a suitable choice of eliminates this undesirable effect.

ii.6 Centering and choice of

The present section essentially relies on formal arguments. To prescribe in Equation (6), we need first to center at the core of the solution of (1) by imposing the supplementary condition


where we have introduced the quantity


and where the constant represents the half-size of the computational box in the numerical calculations. Imposing condition (22) makes the solution of (1) unique (and not only unique up to translations), which is crucial for numerical purposes.

Next, on the basis of expression (17) for , we prescribe the function as


and is a fixed parameter (the reciprocal of some characteristic time). Since by (21) in (6) converges to up to a translation, we formally have by comparing the following expression to (17):


Substituting this limit into definition (24) one deduces that at large times


If we can show that when , which is equivalent to


by definition of , then in the same limit we shall have by (26) . The limit (27) indicates that the choice (24) forces the dynamical solution to obey asymptotically the same centering as . Put differently, this amounts to computing in a comoving frame. The rightmost term in (24) can thus be called a centering correction to the velocity.

To show that , we differentiate the large-time expression (21) of with respect to time. Further invoking (26), we obtain at large times (the dot denotes a total time derivative)


Integrating (28) over then yields the approximate first-order differential equation


where the rightmost expression follows from having neglected the term in the asymptotic expansions (13) of . From (13) this is legitimate if . The latter condition is compatible with being fixed in numerical computations if tends to a finite value at large times. The latter property follows from a simple self-consistent argument. Indeed, Equation (29) implies that for where is some time above which (21) holds, so that effectively . Substituting this expression of into (26), one deduces an approximate analytical expression for that tends to . Writing (21) in the form , and substituting the obtained expression of finally entails the desired saturation property in the form . Unfortunately no quantitative estimate of the latter quantity is available. Yet, on the basis of numerical experiments, it can be expected to be of the order of the (unknown) dislocation core width.

As a final remark, we anticipate by indicating that, later on, will be taken inversely proportional to the algorithmic time step. Thus, in practice in the limit of continuous times. Because of the dependence of , we observe that the centering correction in (24) remains well-behaved in this limit.

ii.7 An analytical solution

For and


with CBCs


the dimensionless Equation (1) admits the following analytical solution, easily deduced from the well-known solution Rosakis () to the original Weertman equation:

with (32b)

This solution is symmetric in the sense of Section II.3. We will use this test case in Sec. V as a benchmark. So, when the core width and the velocity both blow up as


Since and are not the physical velocity and core width, this behavior is not the hallmark of a physical pathology of the model. It however implies that the comptational box should be taken wider and wider to contain the core of , and that the latter moves with nearly infinite velocity, which has numerical consequences to be examined in Sec. V.1.

In the form (13), the asymptotic behaviors deduced from a direct expansion of (32a) read


where the length scales (12) are . Thus, in this particular example where the solution is symmetric the asymptotic behaviors provide a connection between the core width and the next-to-leading terms in the expansion. This supports the interpretation put forward in Sec. II.3 that in generic asymmetric situations represent characteristic scales of variation of the dislocation density.

Iii Building blocks

Our numerical scheme to solve the dynamical system (6) crucially rests on evaluating the action of the operator in the Fourier domain, by means of the Discrete Fourier Transform (DFT). This section explains the underlying spatial discretization procedure, and outlines the key features of the implementation.

iii.1 Temporal and spatial discretization

At discrete times with step , we need a suitable representation of over the whole -axis. To this aim, we define a computational box , discretized into elementary intervals of width . We write the function as


where is a time-evolving part, the support of which is contained within the box , and is a fixed reference function that complies with the asymptotic behaviors (9). The latter function is prescribed once for all, and plays the role of boundary conditions for the operators and . Decomposition (35) is motivated by the fact that the operator is non-local, and that the solution of (1) does not vanish at infinity. Thus, the tail contributions of at infinity should be taken into account when computing . They are represented in (35) by those of . The need to properly account for tail contributions will be illustrated by means of numerical examples in Section V.4.

In this article, we take as the linear combination

where the basis functions are chosen such that can be computed analytically, as

and where is an extra arbitrary scaling parameter. The four coefficients are determined so as to satisfy the four constraints expressed by (9). Comparing the asymptotic expansions of for directly deduced from (36) to the generic expressions (9) leads to


where and and have been defined, respectively, in Equations (14) and (23). Thus, choosing makes vanish. In the exact case of Sec. II.7, where furthermore , only and are nonzero and the solution is already completely retrieved at the level of . Therefore we shall need to take different from to be able to use the exact solution of Sec. II.7 as a non-trivial benchmark of the algorithm in Sec. V. We observe that the asymptotic expansion of in (36) involves only odd inverse powers of .

Finally, introducing discrete positions , the function is represented inside the box by the vector of components . It is decomposed according to (35) as


where the vector of components is now the unknown of the problem.

It must be emphasized that, however convenient, representation (36) is largely arbitrary. Indeed, any other smooth function obeying the asymptotic conditions (9), mostly varying inside the computational box, and such that can be accurately computed once for all (either analytically, or even numerically), would equally well fit our purpose.

iii.2 Numerical Fourier transform and discretization of the diffusion operator

This section addresses the discretization of the linear integro-differential operator . In view of (2), the latter involves a convolution by a pseudofunction KANW04 (), which can be done in the Fourier representation. Crucially, the present approach uses the continuous Fourier form of the operator because this representation is straightforwardly diagonal, and versatile in the sense that it could as well be employed to address the fractional Laplacian .

The FT will be implemented in DFT form, which allows one to benefit from Fast-Fourier-Transform (FFT) routines. As is well known, carrying out a convolution over the real axis by means of DFT and multiplications in the Fourier domain turns this convolution into a periodic one, which induces undesirable periodicity effects in the direct space near both extremities of the interval of interest Bracewell ().

The usual workaround is to use zero-padding. Briefly, the technique consists in doubling the spatial extent of the interval of interest, in such a way that the undesirable effects that take place during convolution remain confined to the extra region. The vector to be convolved is continued by zero in the latter region. Once the convolution with the kernel has been carried out in the Fourier space by Fourier component multiplication, and the result has subsequently been transformed back to direct space by DFT inversion, the irrelevant contributions of the added region that have been spoiled by periodicity effects are dropped. This last step constitutes a projection of the result onto the initial interval of interest.

For details the reader is referred to the classical reference (NumRecipes, , p. 643). However, the latter reference only addresses situations where the initial vector and the convolution kernel are both known in the direct space. In contrast, since the present approach uses the convolution kernel in continuous Fourier form, it requires identifying the Fourier modes in the continuum to the discrete modes used by the DFT. For that reason, the labelling and ordering of the vector components hereafter somewhat differ from the ones in Ref. NumRecipes (). However, as Equation (46) below shows, the final outcome will take the form of a standard convolution in the direct space.

Thus, zero-padding is carried out by means of the injection that maps into , and the inverse transform involves a projector that maps onto . These operators are defined by


Note that leaves invariant.

The DFT that operates on extended vectors is denoted hereafter by (or, alternatively, by ). The corresponding enlarged spatial grid that extends the one in Equation (38) is for . It admits the dual wavemode grid


Accordingly, DFTs are carried out for each wave mode as


From (35) we can now compute the discretized form of as


in which is evaluated at point by means of (36), and is computed via the DFT as


Combining (43) and (44) we shall introduce , the discretization of , defined as


A straightforward check shows that the rightmost term simplifies as


Although there are wave modes, the sum only involves indices , precisely because of the use of injection . As desired, the convolution in (46) is well-behaved and not spoiled by undesirable periodic effects, since it distinguishes at each point the contributions of the on the left of , from those on the right for . This would not have been the case, had fewer points been retained to compute DFTs.

iii.3 Discretization of the advection operator

We discretize the advection operator with the following upwind scheme Fletcher ():


where and , and where the operators are implemented via expressions with error (see (Fletcher, , p. 297) with ), namely,


The choice a third-order implementation is motivated in Sec. V.4 by numerical considerations (see Ref. Fletcher () for schemes of orders and and Ref. (Hildebrand, , p. 111) for order ). Since


the discretization (47) of the advection operator actually reads


where only depends on the sign of . In this expression, plays the role of boundary conditions when computing at the extremities of the computational box, e.g., for in (50).

Note that whereas, in principle, implementing should require adding only a few extra points at the exterior of the interval , our implementation (50) is carried out in practice on points, due to the introduction of and for consistency with the implementation of the operator . Also, remark that the operators are almost diagonal, in the sense that their periodizations are diagonal in the Fourier space. In this connection, we introduce for further use in Section IV.2 the vectors defined as


with , and where the periodic convention applies in the evaluation of the operator in (51).

iii.4 Velocity computation

As seen above, Equation (6) can be solved in the comoving frame by using the velocity given by (24). In this way, the dislocation core lies as remote as possible from the box boundaries to minimize the influence of the approximations made in handling the tails. To proceed, we substitute in (50), at each time , the quantity by (with the convention that ) computed from a discretized version over points of expression (24), in which has been replaced by (since and ); namely,


In these expressions , and are the weights , , , …, , , of the Simpson integration rule. Remark that depends only on via its sign, which is of little consequence except in calculations at small where is close to and may oscillate during iterations. The term within brackets results from a straightforward discretization of the integral in (24), in which the tail contributions have been evaluated analytically from the asymptotic expansions of (for simplicity, we have not evaluated the full tail contributions of ).

Moreover, a suitable value of stems from the empirical consideration that if we discretize Equation (29) in explicit Euler form as , monotone convergence towards 0 of is ensured by taking . Correspondingly, the value of used henceforth in (52a) is .

Iv Algorithm

This section describes the iterative numerical scheme used to compute and . This algorithm is explicit in time, is consistent with the dynamical system (6), combines the above-discretized operators, and remains stable even with a large time step .

iv.1 Procedure

Computations go as follows. First, for given local minimizers and of , the algorithm is typically initialized by choosing arbitrarily the values inside the box, preferentially not too far from the expected solution. This can be done in a number of ways; notably, by using as an initial condition the function with chosen large enough to encompass the expected overall width of the dislocation density (typically, a few times the characteristic scale ), whence . A few low-resolution runs may help adjusting . Obviously, to get a reasonable representation of the solution, the discretization step must necessarily be less than . Also, when carrying out incremental parametric studies, the solution computed from the previous value of the parameter under consideration can be used as an initial condition, to save CPU time. However, for studying stability issues, we shall purposely take initial data far from the expected shape of a dislocation.

Denoting by the scheme presented below, we iterate


until the difference between the results of two successive iterations is small, in the sense that


where is a user-defined stopping criterion. Upon completion at some , the algorithm returns and the associated , evaluated thanks to (52a). Unless otherwise stated, in the numerical calculations of Sec. V below.

iv.2 The Preconditioned Collocation Scheme

The scheme described hereafter, which we call the Preconditioned Collocation Scheme (PCS), is based on the requirement that the long-time limit of should solve the following static equation:


which is the discretized form of (1). This is a collocation method. A first naive way to proceed would be to attempt solving (6) by writing


where should be adjusted in order to achieve convergence. At convergence, the solution obeys (55) and depends on only through the evaluation of the numerical velocity defined by (52a). As will be justified in Section V.3 this dependence is of little relevance, which is an advantage of this approach.

However, system (56) is ill-conditioned, because it involves the stiff operators and