0905/0905.4617.md

Phase-Field Modeling of Nonlinear Material Behavior

Y.-P. Pellegrini, C. Denoual and L. Truskinovsky

Abstract Materials that undergo internal transformations are usually described in solid mechanics by multi-well energy functions that account for both elastic and transformational behavior. In order to separate the two effects, physicists use instead phase-field-type theories where conventional linear elastic strain is quadratically coupled to an additional field that describes the evolution of the reference state and solely accounts for nonlinearity. In this paper we propose a systematic method allowing one to split the nonconvex energy into harmonic and nonharmonic parts and to convert a nonconvex mechanical problem into a partially linearized phase-field problem. The main ideas are illustrated using the simplest framework of the Peierls–Nabarro dislocation model.

1 Introduction

Nonconvex energy potentials are used in solid mechanics for the modeling of martensitic transformations [9], plasticity [1] and fracture [25]. Parts of the resulting energy landscapes correspond to sufficiently smooth deformations preserving the locally affine structure of the lattice environment of each atom. Other parts represent highly distorted atomic arrangements associated with either loss or reacquisition of nearest neighbors. While deformations of the first type can (often) be described by

---

Y.-P. Pellegrini and C. Denoual

CEA, DAM, DIF, F-91297 Arpajon, France; e-mail: yves-patrick.pellegrini@cea.fr,
christophe.denoual@cea.fr

L. Truskinovsky

Laboratoire de Mécanique des Solides, CNRS UMR-7649, École Polytechnique, Route de Saclay,
F-91128 Palaiseau Cedex, France; e-mail: trusk@lms.polytechnique.fr

In: K. Hackl (ed.) Proceedings of the IUTAM Symposium on Variational Concepts with Applications to the Mechanics of Materials, September 22–26, 2008, Bochum, Germany (Springer-Verlag, presumably 2010). In press.the conventional strain tensor of (linear) elasticity theory, a representation of the deformations of the second type requires introducing additional internal variables accounting for deviations from the local affinity of the stressed atomic configurations. In particular, these supplementary variables describe the evolution of the local reference state (LRS) from which the elastic deformations are measured [3, 11, 26]. The main difference between the elastic strains and these supplementary internal variables is that the dynamics of the former is typically inertial, while that of the latter is usually overdamped. Sometimes the nonelastic variables can be minimized out as in the case of deformational plasticity (e.g., [2]). In this paper we deal instead with situations where the internal variables have to be revealed rather than hidden.

We assume that the coarse-grained nonconvex energy density $f(\varepsilon)$ is known either from extrapolations of experimental measurements or from ab-initio calculations involving atomic homogeneity constraints. We suppose that the argument $\varepsilon$ of this function, that represents a coarse-grained strain, is small and can be additively split into the linear elastic part $e$ , and a phase-field part $\eta$ that accounts for the nonelastic evolution of the LRS. Our next assumption is that $f$ can be represented as a sum of two terms: the elastic energy $f_e$ , which depends on $e = \varepsilon - \eta$ and the phase-field energy $g$ , which depends on $\eta$ . We interpret $f(\varepsilon)$ as the outcome of adiabatic elimination of the variable $\eta$ and consider the inverse problem of recovering the phase-field energy $g(\eta)$ from the function $f(\varepsilon)$ under the assumption that the function $f_e(e)$ is quadratic. The problem of the identification of $g(\eta)$ reduces to a problem of optimization and the relation between the 'optimally' related functions $f(\varepsilon)$ and $g(\eta)$ is studied in some prototypical cases. If, in contrast, the function $g(\eta)$ is chosen independently, the corresponding function $f(\varepsilon)$ is typically non-smooth and non single-valued, e.g. [6].

To motivate the need for the phase-field variables we consider in full detail a specific physical example. It deals with the mixed, discrete-continuum representation of a dislocation core [12, 16]. More specifically, we develop a modified version of the classical Peierls–Nabarro (PN) model that accounts for a finite thickness of the slip region. In this problem the coarse-grained description of the slip zone is provided by the so-called $\gamma$ -potential [5, 27]. The phase field represents an “atomically sharp” slip and the part of the interaction potential related to $g$ gives rise to the slip-related pull-back force [7, 16, 23]. Our general method of recovering the expression for this force represents an extension of Rice’s transform, which was first introduced in the context of a dislocation nucleation problem [19].

In this paper only the simplest scalar problem in a one-dimensional setting is considered. The slightly more general question of extracting from the coarse-grained energy a convex (instead of quadratic) component will be examined elsewhere [17].

2 Surface Problem

We begin with the special case when the phase field is localized on a surface. In problems involving fracture or slip it often proves convenient to represent the energyof a body as the sum of a bulk term depending on strain gradients and a surface term penalizing displacement discontinuities. The bulk term is usually modeled by linear elasticity. The modeling of the surface energy is less straightforward [6, 25]. For instance, the models will be different depending on whether the location of the discontinuities is known a priori or not.

In a 1D setting with a known fracture set the equilibrium problem reduces to minimizing the following energy functional

$$W[u]={\int }_0^{1}dx{f}_e(u_x)+\sum _{{\mathrm{\Gamma }}_a}{f}_a(\delta (x))\mathrm{.}(1)W[u]\; =\; \backslash int\_0^1\; dx\; f\_e(u\_x)\; +\; \backslash sum\_\{\backslash Gamma\_a\}\; f\_a(\backslash delta(x)).\; \backslash quad\; (1)$$

Here $f_e(e) = (E/2)e^2$ , $E > 0$ is the elastic modulus and $a$ is a coarse-graining length scale that typically exceeds several atomic sizes. The set $\Gamma_a$ in (1) represents discontinuity points resolved at scale $a$ and $\delta(x) = [u]_a(x)$ is the corresponding displacement discontinuity. The surface energy $f_a(\delta)$ is then an effective interaction over the distance $a$ ; in particular, the shear-related component of $f_a(\delta)$ coincides the $\gamma$ -potential mentioned in the Introduction.

In the case when the fracture set is unknown the surface energy has to be chosen differently. The reason is that in this model the displacement discontinuity at scale $a$ does not represent the microscopic slip between neighboring atomic planes, and therefore the difference between elastic deformation and inelastic slip has to be yet resolved at this scale [19]. More precisely, linear elasticity, which has nothing to do with slip and which is already accounted for in the bulk term, has not been excluded from $f_a(\delta)$ . The identification of the surface energy in (1) with $f_a(\delta)$ , which is quadratic at the origin, leads in a free discontinuity problem to a degenerate solution with infinitely many infinitely small discontinuities [6].

To remove linear elasticity from the surface term, one should replace the coarse-grained discontinuity $[u]_a$ by the atomically sharp slip $\eta(x) = u$ that does not depend on $a$ . The energy (1) is then rewritten as

$$W[u]={\int }_0^{1}dx{f}_e(u_x)+\sum _{\mathrm{\Gamma }}g(\eta ),(2)W[u]\; =\; \backslash int\_0^1\; dx\; f\_e(u\_x)\; +\; \backslash sum\_\{\backslash Gamma\}\; g(\backslash eta),\; \backslash quad\; (2)$$

where now $\Gamma$ is the set of discontinuity points corresponding to $a = 0$ . The problem is to find the relation between the function $f_a(\delta)$ , representing an empirical input, and the unknown function $g(\eta)$ .

To define $g(\eta)$ we divide the total slip $\delta$ into an elastic part, $ae$ , where $e$ is an equivalent elastic strain, and an inelastic part $\eta$ . The function $g(\eta)$ is defined by the condition that $f_a(\delta)$ is a relaxation of the energy $af_e(e) + g(\eta)$ under the condition that $ae + \eta = \delta$ , namely:

$$f_a(\delta )=\underset{\eta }{\inf }[a\frac{E}{2}{\left(\frac{\delta -\eta }{a}\right)}^2+g(\eta )]\mathrm{.}(3)f\_a(\backslash delta)\; =\; \backslash inf\_\{\backslash eta\}\; \backslash left[\; a\; \backslash frac\{E\}\{2\}\; \backslash left(\; \backslash frac\{\backslash delta\; -\; \backslash eta\}\{a\}\; \backslash right)^2\; +\; g(\backslash eta)\; \backslash right].\; \backslash quad\; (3)$$If the energy $f_a(\delta)$ is a single-well function and the infimum is unique, the function $g$ is completely defined. If $f_a(\delta)$ is periodic as in the case of dislocations, in order to have a uniquely-defined $g(\eta)$ , we need to replace in definition (3) the global minimization by a properly-defined local minimization denoted hereafter by 'inf_loc' (minimization over $\eta$ starting from the minimum of $f_a$ closest to $\delta$ .)

In what follows, our task will be to reverse definition (3) and to recover the nonequilibrium energy $g(\eta)$ from $f_a(\delta)$ . What allows us to proceed is the specific (harmonic) structure of the elastic part of the energy.

We observe that the function $g$ must satisfy the following necessary condition

$$(E\mathrm{/}a)(\delta -\eta )=g^{\mathrm{\prime }}(\eta )\mathrm{.}(4)(E/a)(\backslash delta\; -\; \backslash eta)\; =\; g\text{'}(\backslash eta).\; \backslash quad\; (4)$$

Moreover, differentiation of (3) wrt. $\delta$ gives

$$f_a^{\mathrm{\prime }}(\delta )=(E\mathrm{/}a)(\delta -\eta )\mathrm{.}(5)f\text{'}\_a(\backslash delta)\; =\; (E/a)(\backslash delta\; -\; \backslash eta).\; \backslash quad\; (5)$$

These two equations allow one to represent $g'(\eta)$ in the following parametric form [7, 8, 19]

$$(\eta ,g^{\mathrm{\prime }}(\eta ))=(\delta -(a\mathrm{/}E)f_a^{\mathrm{\prime }}(\delta ),f_a^{\mathrm{\prime }}(\delta ))\mathrm{.}(6)(\backslash eta,\; g\text{'}(\backslash eta))\; =\; (\backslash delta\; -\; (a/E)f\text{'}\_a(\backslash delta),\; f\text{'}\_a(\backslash delta)).\; \backslash quad\; (6)$$

The parametric representation for $g(\eta)$ then reads

$$(\eta ,g(\eta ))=(\delta -(a\mathrm{/}E)f_a^{\mathrm{\prime }}(\delta ),f_a(\delta )-[a\mathrm{/}(2E)][f_a^{\mathrm{\prime }}(\delta ){]}^2)\mathrm{.}(7)(\backslash eta,\; g(\backslash eta))\; =\; \backslash left(\; \backslash delta\; -\; (a/E)f\text{'}\_a(\backslash delta),\; f\_a(\backslash delta)\; -\; [a/(2E)]\; [f\text{'}\_a(\backslash delta)]^2\; \backslash right).\; \backslash quad\; (7)$$

Since for nonconvex $f_a(\delta)$ this representation may lead to a multivalued function $g(\eta)$ formula (7) must be supplemented by an additional branch selection procedure. To illustrate the mapping $f(\delta) \rightarrow g(\eta)$ given by (7) and the selection of a

Fig. 1 Parametric transforms (6), (7) applied to a Lennard-Jones potential $f(\delta) = \delta^{-12} - \delta^{-6}$ . a) $f(\delta)$ (dashed), and the two branches of $g(\eta)$ (solid), where 'p' ('u') labels the physical (unphysical) branch; b) $f'(\delta)$ and the 'p' branch of $g'(\eta)$ ; c) $\eta(\delta)$ .

physical branch we consider a Lennard-Jones potential $f_a$ , with $a = 1$ and assume that $E = f''(\delta_0)$ , where $\delta_0$ is the only minimum of $f$ (see Fig. 1). Notice that theresulting function $g'(\eta)$ has an infinite slope at $\eta = \delta_0$ and that for $\eta \gtrsim \delta_0$ we must have $g(\eta) \propto (\eta - \delta_0)^{3/2}$ .

The removal of the linear elastic part of the energy becomes important in PN-type modeling of dislocations. Consider, for instance, a straight screw dislocation in an isotropic linear-elastic body and assume that the sharp discontinuity plane, $y = 0$ , lies between the two effective gliding surfaces located at $y = \pm a/2$ . To account for the finite thickness of the core region $a$ we need to modify the classical PN model [12]. According to our interpretation the linear elastic stress outside the slip region $(-a/2, a/2)$ must be balanced by the coarse-grained pull-back stress that is resolved at the spatial scale $a$ . We therefore interpret the pull-back stress at this scale as $f'_a(\delta(x))$ , where $f_a$ is the $\gamma$ -potential, a periodic function with period $b$ and with $f'_a(0) = 0$ . The expression for the linear stress outside the slip region is derived in the Appendix. With these considerations in mind we obtain for the unknown function $\eta(x)$ representing a mathematical slip at $y = 0$ the following system of equations

$$-\frac{\mu }{\pi a}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}d{x}^{\mathrm{\prime }}{\eta }^{\mathrm{\prime }}(x^{\mathrm{\prime }})\arctan \frac{a}{2(x-x^{\mathrm{\prime }})}+{\bar{\sigma }}_a(x)=f_a^{\mathrm{\prime }}(\delta (x)),(8)-\backslash frac\{\backslash mu\}\{\backslash pi\; a\}\; \backslash int\_\{-\backslash infty\}^\{+\backslash infty\}\; dx\text{'}\; \backslash eta\text{'}(x\text{'})\; \backslash arctan\; \backslash frac\{a\}\{2(x-x\text{'})\}\; +\; \backslash bar\{\backslash sigma\}\_a(x)\; =\; f\text{'}\_a(\backslash delta(x)),\; \backslash quad\; (8)$$

$$\delta (x)=(a\mathrm{/}\mu )f_a^{\mathrm{\prime }}(\delta (x))+\eta (x),(9)\backslash delta(x)\; =\; (a/\backslash mu)\; f\text{'}\_a(\backslash delta(x))\; +\; \backslash eta(x),\; \backslash quad\; (9)$$

where $\bar{\sigma}_a$ is the resolved applied stress at scale $a$ . If we match the linear elastic behavior at $\eta = 0$ with that in the bulk regions we obtain that $\mu = a f''_a(0)$ . Using in this relation the physical shear modulus and the value of $f''_a(0)$ from the $\gamma$ -potential provides a rough estimate for $a$ , the effective interaction range.

We notice that parameter $a$ enters both equations (8,9), which makes this system different from the one studied in [16, 19]. The ideas behind our nonlocal extension of the PN model are also different from that of Ref. [13] where a nonlocal kernel was introduced empirically as part of the pull-back stress, and the usual $1/(x-x')$ kernel was used for the bulk stress.

To bring the system (8,9) into the framework of phase-field models, we identify the effective pull-back force $f'_a(\delta(\eta))$ with $g'(\eta)$ and rewrite Eq. (8) as

$$\frac{\mu }{2}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}d{x}^{\mathrm{\prime }}{K}_a(x-x^{\mathrm{\prime }}){\eta }^{\mathrm{\prime }}(x^{\mathrm{\prime }})+{\bar{\sigma }}_a(x)=g^{\mathrm{\prime }}(\eta )\mathrm{.}(10)\backslash frac\{\backslash mu\}\{2\}\; \backslash int\_\{-\backslash infty\}^\{+\backslash infty\}\; dx\text{'}\; K\_a(x-x\text{'})\; \backslash eta\text{'}(x\text{'})\; +\; \backslash bar\{\backslash sigma\}\_a(x)\; =\; g\text{'}(\backslash eta).\; \backslash quad\; (10)$$

where $K_a(x) = -(2/\pi a) \arctan(a/2x)$ . It is now easy to see that $g'(\eta)$ enjoys the parametric representation

$$(\eta ,g^{\mathrm{\prime }}(\eta ))=(\delta -\frac{a}{\mu }{f}_a^{\mathrm{\prime }}(\delta ),f_a^{\mathrm{\prime }}(\delta )),(11)(\backslash eta,\; g\text{'}(\backslash eta))\; =\; \backslash left(\; \backslash delta\; -\; \backslash frac\{a\}\{\backslash mu\}\; f\text{'}\_a(\backslash delta),\; f\text{'}\_a(\backslash delta)\; \backslash right),\; \backslash quad\; (11)$$

where we recognize the mapping (6) (see also [16, 19, 23]). To make the link with the classical PN model one needs to consider the limit $a \rightarrow 0$ . By computing $\eta$ in terms of $\delta$ and expanding (10) in powers of $a$ , we obtain to order $O(a)$ the following ‘gradient’ extension of the PN model$$-\frac{\mu}{2\pi} \int_{-\infty}^{+\infty} dx' \frac{\delta'(x')}{x-x'} + \lambda \delta''(x) + \tilde{\sigma}_a(x) = f(\delta(x)), \quad (12)$$

where $\lambda = a\mu/4$ . For different weakly or strongly nonlocal generalizations of the PN model see [13, 21]. Equation (12) features an effective applied stress that differs from $\sigma_a(x, 0)$ defined in the Appendix by an $O(a)$ correction, namely $\tilde{\sigma}_a(x) \equiv \sigma_a(x, 0) + (a/2)[(1/2)\partial_y \sigma_a(x, 0) - K_0 \star \partial_x \sigma_a(x, 0)]$ . The classical PN model is retrieved by letting $a = 0$ .

3 Bulk Problem

Now let us place the problem in a more general framework. The task is to approximate locally the empirical potential $f(\varepsilon)$ by a quadratic function with an optimally chosen reference state $\eta$ , and to associate with this state a reference energy $g(\eta)$ . Behind such construction is the assumption that all the nonlinearity of the problem is related to the evolution of the reference state. The simplest setting to pose formally the problem is the one-dimensional geometrically linearized theory of nonlinear elastic bars.

According to our interpretation the empirical energy is represented as

$$f(\epsilon )=\underset{\eta ,\text{loc}}{\inf }[\frac{E}{2}(\epsilon -\eta {)}^2+g(\eta )](13)f(\backslash varepsilon)\; =\; \backslash inf\_\{\backslash eta,\; \backslash text\{loc\}\}\; \backslash left[\; \backslash frac\{E\}\{2\}\; (\backslash varepsilon\; -\; \backslash eta)^2\; +\; g(\backslash eta)\; \backslash right]\; \backslash quad\; (13)$$

and the problem is to find the intrinsic phase-field function $g(\eta)$ . Following the previous section we write the parametric representation for $g'(\eta)$ in the form

$$(\eta ,g^{\mathrm{\prime }}(\eta ))=(\epsilon -\frac{f^{\mathrm{\prime }}(\epsilon )}{E},f^{\mathrm{\prime }}(\epsilon ))\mathrm{.}(14)(\backslash eta,\; g\text{'}(\backslash eta))\; =\; \backslash left(\; \backslash varepsilon\; -\; \backslash frac\{f\text{'}(\backslash varepsilon)\}\{E\},\; f\text{'}(\backslash varepsilon)\; \backslash right).\; \backslash quad\; (14)$$

The function $g(\eta)$ is then given by the mapping

$$(\eta ,g(\eta ))=(\epsilon -\frac{f^{\mathrm{\prime }}(\epsilon )}{E},f(\epsilon )-\frac{f^{\mathrm{\prime }}(\epsilon {)}^2}{2E})\mathrm{.}(15)(\backslash eta,\; g(\backslash eta))\; =\; \backslash left(\; \backslash varepsilon\; -\; \backslash frac\{f\text{'}(\backslash varepsilon)\}\{E\},\; f(\backslash varepsilon)\; -\; \backslash frac\{f\text{'}(\backslash varepsilon)^2\}\{2E\}\; \backslash right).\; \backslash quad\; (15)$$

The consistency of this procedure requires the parameter $E$ and the function $f(\varepsilon)$ to be related. If we expand the parametric definition of $g(\eta)$ near a reference state $\varepsilon_0$ where $f'(\varepsilon_0) = 0$ , we obtain $g(\varepsilon_0) = f(\varepsilon_0)$ , $g'(\varepsilon_0) = 0$ and $g''(\varepsilon_0) = f''(\varepsilon_0)/[1 - f''(\varepsilon_0)/E]$ . The natural choice $E = f''(\varepsilon_0)$ makes $g''(\varepsilon_0)$ infinite. The behavior of the higher derivatives of $g(\eta)$ near $\eta = \varepsilon_0$ is found by assuming (without loss of generality) that derivatives $f^{(k)}(\varepsilon_0)$ vanish for $k = 3, \dots, n-1$ . The order of the asymptotics depends on $n > 2$ , which is the first integer such that $f^{(n)}(\varepsilon_0) \neq 0$ :

$$(\eta ,g(\eta ))\simeq ({\epsilon }_0-\frac{f^{(n)}({\epsilon }_0)}{(n-1)!E}\delta {\epsilon }^{(n-1)},f({\epsilon }_0)-\frac{(n-1)f^{(n)}({\epsilon }_0)}{n!}\delta {\epsilon }^n)\mathrm{.}(16)(\backslash eta,\; g(\backslash eta))\; \backslash simeq\; \backslash left(\; \backslash varepsilon\_0\; -\; \backslash frac\{f^\{(n)\}(\backslash varepsilon\_0)\}\{(n-1)!E\}\; \backslash delta\; \backslash varepsilon^\{(n-1)\},\; f(\backslash varepsilon\_0)\; -\; \backslash frac\{(n-1)f^\{(n)\}(\backslash varepsilon\_0)\}\{n!\}\; \backslash delta\; \backslash varepsilon^n\; \backslash right).\; \backslash quad\; (16)$$Hence $g$ behaves near its minimum as: $|g(\eta) - g(\varepsilon_0)| \sim |\eta - \varepsilon_0|^{n/(n-1)}$ . The generic case is $n = 3$ ; the case $n = 4$ corresponds to periodic potential relevant for dislocations; for $f$ locally harmonic, $n = \infty$ .The working of Eqs. (14-17) with $E = f''(\varepsilon_0)$ is illustrated in Fig. 2. In the domain at the left of $\varepsilon_0$ , where $f$ grows faster than harmonic the desired tangency point does not exist. In this case the difference $f(\varepsilon) - \frac{E}{2}(\varepsilon - \eta)^2$ is maximized at $\varepsilon = -\infty$ . This situation takes place in the Lennard-Jones example of Sec. 2 where we have to use $g(\delta) = +\infty$ for $\delta < \delta_0$ (hatched area of Fig. 1a).

To handle general multi-well energies, we first introduce the Stillinger-Weber mapping $\varepsilon_0(\varepsilon)$ that links to any state $\varepsilon$ the local minimum $\varepsilon_0$ of $f(\varepsilon)$ that would be attained from this state by steepest-descent [22]. Next, we modify equations (13) and (17) as:

$$f(\epsilon )=\underset{\eta ,\text{loc}}{\inf }[\frac{1}{2}{f}^{\mathrm{\prime }\mathrm{\prime }}({\epsilon }_0(\epsilon ))(\epsilon -\eta {)}^2+g(\eta )],(18)f(\backslash varepsilon)\; =\; \backslash inf\_\{\backslash eta,\; \backslash text\{loc\}\}\; \backslash left[\; \backslash frac\{1\}\{2\}\; f\text{'}\text{'}(\backslash varepsilon\_0(\backslash varepsilon))\; (\backslash varepsilon\; -\; \backslash eta)^2\; +\; g(\backslash eta)\; \backslash right],\; \backslash quad\; (18)$$

$$g(\eta )=\underset{\epsilon ,\text{loc}}{\sup }[f(\epsilon )-\frac{1}{2}{f}^{\mathrm{\prime }\mathrm{\prime }}({\epsilon }_0(\eta ))(\epsilon -\eta {)}^2]\mathrm{.}(19)g(\backslash eta)\; =\; \backslash sup\_\{\backslash varepsilon,\; \backslash text\{loc\}\}\; \backslash left[\; f(\backslash varepsilon)\; -\; \backslash frac\{1\}\{2\}\; f\text{'}\text{'}(\backslash varepsilon\_0(\backslash eta))\; (\backslash varepsilon\; -\; \backslash eta)^2\; \backslash right].\; \backslash quad\; (19)$$

Whereas (19) defines $g$ , equation (18) states that knowing $f$ is equivalent to knowing $g$ plus the linear-elastic behavior of $f$ near its local minima.

The precise meaning of the “loc” in Eqs. (18), (19) is as follows. Operationally, the minimization in the definition of $f$ is carried out over $\eta$ , starting from $\varepsilon_0(\varepsilon)$ , the local minimum nearest to $\varepsilon$ determined by the SW mapping; the corresponding elastic modulus is also determined by the starting point. The maximization in the definition of $g(\eta)$ proceeds along similar lines except that now the relevant elastic modulus is determined by the local minimum closest to $\eta$ , and is fixed during the maximization. The optimization is carried out starting from $\varepsilon = \eta$ .

Figs. 3 illustrate the case of a double-well potential with unequal curvatures of the wells. Notice that in contrast to what we saw in Fig. 1 the function $\eta(\varepsilon)$ is now bounded.

Fig. 4 Transform (19) applied to the piecewise-harmonic periodic potential $f(\delta) = (\delta - [\delta])^2/2$ . Notation $[\cdot]$ stands for the integer part. The function $\eta(\delta)$ is determined as the argmin in (18).

Another interesting case is the periodic potential that is used in the description of reconstructive phase transitions (e.g., [4]). Consider, for instance, the piecewise-harmonic periodic case shown in Fig. 4 that is often used in analytical studies [13, 21]. The parametric representation (14) of $g$ is here useless and the definition (19) must be used instead. In this extreme case, all elasticity has been removed from $g$ and the resulting $g(\eta)$ is cone-shaped at its minima (Fig. 4a) as predicted by Eq. (16) for $n \rightarrow \infty$ . The force $g(\eta)$ is discontinuous (Fig. 4b) and its extreme values provide thresholds for the evolution of $\eta$ , whose stepwise character is an artifact due to the absence of smooth spinodal regions in $f$ .

It is also instructive to consider for comparison the case of an unbounded harmonic potential $f(\varepsilon) = (E_f/2)(\varepsilon - \varepsilon_0)^2$ . From (17) with $E = E_f$ , one deduces that $g(\eta) = +\infty$ if $\eta \neq \varepsilon_0$ and $g(\eta = \varepsilon_0) = 0$ . This trivial example indicates that in a purely linear-elastic model, the reference state does not have to evolve.

4 Concluding Remarks

The goal of this paper was to reveal in the simplest setting the variational nature of the generalized Rice transform. The problem consists in splitting a coarse-grained lattice potential $f$ , describing the overall deformation of a sufficiently large number of atoms, into a (quasi) convex elastic potential and an inelastic potential $g$ dealing with structural rearrangements. Here the potential $f$ is assumed to be measurable by molecular statics along a prescribed deformation path relevant to the material transformation in question. For simplicity, the elastic potential is assumed in this paper to be a standard quadratic function of the macroscopic strain. The inelastic potential $g$ must be a function of the phase-field variable $\eta$ , whose identification represents an important part of the problem. While our precise construction solving the above problem is presented in the static setting (see Eq. (19)), the motivation for the splitting concerns, first of all, dynamical applications (e.g., [7]). Thus we assume that material displacement $u$ associated to the strain $\varepsilon = \partial_x u$ evolves inertially almost without damping (standard elastodynamics), while the dynamics of the phase-field variable is overdamped. More precisely we assume that the relaxation of the variable $\eta$ follows the time-dependent Ginzburg-Landau (TDGL) equation. By means of an empirical ‘viscosity’ parameter $\nu$ we can write the evolution equation in the form

$$\dot{\eta }=-\frac{1}{\nu }\frac{\mathrm{\partial }}{\mathrm{\partial }\eta }[\frac{1}{2}{f}^{\mathrm{\prime }\mathrm{\prime }}({\epsilon }_0(\epsilon ))(\epsilon -\eta {)}^2+g(\eta )],\backslash dot\{\backslash eta\}\; =\; -\backslash frac\{1\}\{\backslash nu\}\; \backslash frac\{\backslash partial\}\{\backslash partial\; \backslash eta\}\; \backslash left[\; \backslash frac\{1\}\{2\}\; f\text{'}\text{'}(\backslash varepsilon\_0(\backslash varepsilon))\; (\backslash varepsilon\; -\; \backslash eta)^2\; +\; g(\backslash eta)\; \backslash right],$$

where we have omitted for simplicity the conventional gradient-penalizing terms (e.g., [7, 24]). In the static setting the above equation reduces to our basic Equ. (18). The definitive advantage of separating the wave motion from an overdamped TDGL relaxation is the possibility to attribute effective damping only to large atomic displacements. Our preliminary investigations [17] indicate that extending the variational set-up presented in this paper to higher dimensions and generalizing it in the direction of extracting (quasi) convex, rather than merely quadratic elastic compo-nents, is feasible. These issues will be addressed systematically in a separate publication.

Appendix

The following computations are largely based on the Eshelby's arguments presented in [10]. Consider a Volterra screw dislocation with zero-width core and with Burgers vector $b$ . The displacement $u_z(x, y)$ has the form

$$u_z(x,y)=\frac{b}{2\pi }\text{Arg}(x+iy)=\frac{b}{2\pi }\arctan \frac{y}{x}+\frac{b}{2}\text{sign}(y)\theta (-x),(20)u\_z(x,\; y)\; =\; \backslash frac\{b\}\{2\backslash pi\}\; \backslash text\{Arg\}(x\; +\; iy)\; =\; \backslash frac\{b\}\{2\backslash pi\}\; \backslash arctan\; \backslash frac\{y\}\{x\}\; +\; \backslash frac\{b\}\{2\}\; \backslash text\{sign\}(y)\; \backslash theta(-x),\; \backslash quad\; (20)$$

where $\theta$ is the Heaviside function, and where the indeterminacy in the discontinuity of $u_z$ is resolved by specifying the glide plane ( $y = 0$ ). The distributional part in the r.h.s. of Eq. (20), usually omitted in the literature (e.g. [12]), is crucial to the present derivation because it represents the irreversible atomic displacements on the plane $y = 0$ . We introduce the eigendistortion, $\beta_{ij}^*$ , as the part of the dislocation-induced distortion $\beta_{ij} = u_{i,j}$ , that is not linear-elastic. For our dislocation, its only non-zero component is $\beta_{yz}^*(x, y) = b \theta(-x) \delta_D(y)$ , where $\delta_D$ is the Dirac distribution [14]. The linear-elastic distortion, $\beta_{ij}^e$ , is defined through the additive decomposition of the total distortion $\beta_{ij}$ , namely $\beta_{ij}^e \equiv \beta_{ij} - \beta_{ij}^*$ [14]. The elastic strains are $e_{ij} = \text{sym} \beta_{ij}^e$ . By using the identity $[\arctan(1/x)]' = \pi \delta_D(x) - 1/(1+x^2)$ , we obtain that the distributional parts in $\beta_{ij}$ and $\beta_{ij}^*$ mutually cancel out giving the standard result [12]

$$e_{xz}(x,y)=-\frac{b}{4\pi }\frac{y}{x^2+y^2},e_{yz}(x,y)=\frac{b}{4\pi }\frac{x}{x^2+y^2}\mathrm{.}(21)e\_\{xz\}(x,\; y)\; =\; -\backslash frac\{b\}\{4\backslash pi\}\; \backslash frac\{y\}\{x^2\; +\; y^2\},\; \backslash quad\; e\_\{yz\}(x,\; y)\; =\; \backslash frac\{b\}\{4\backslash pi\}\; \backslash frac\{x\}\{x^2\; +\; y^2\}.\; \backslash quad\; (21)$$

The stress induced by the eigenstrain is then $\sigma_{iz}^*(x, y) = 2\mu e_{iz}(x, y)$ , where $i = x, y$ .

In the presence of an applied shear stress [15] $\sigma_a \equiv \sigma_{ayz}$ , Eq. (20) becomes

$$u_z(x,y)=\frac{1}{\mu }{\int }_0^{y}d{y}^{\mathrm{\prime }}{\sigma }_a(x,y^{\mathrm{\prime }})+\frac{b}{2\pi }\arctan \frac{y}{x}+\frac{b}{2}\text{sign}(y)\theta (-x),(22)u\_z(x,\; y)\; =\; \backslash frac\{1\}\{\backslash mu\}\; \backslash int\_0^y\; dy\text{'}\; \backslash sigma\_a(x,\; y\text{'})\; +\; \backslash frac\{b\}\{2\backslash pi\}\; \backslash arctan\; \backslash frac\{y\}\{x\}\; +\; \backslash frac\{b\}\{2\}\; \backslash text\{sign\}(y)\; \backslash theta(-x),\; \backslash quad\; (22)$$

The total stress is then $\sigma = \sigma^* + \sigma_a$ and $e = \sigma/2\mu$ . Now, the key step consists in averaging the stress over the layer of width $a$ containing the glide plane. Introduce: $\bar{\sigma}{ij}(x) \equiv \frac{1}{a} \int{-a/2}^{+a/2} dy \sigma_{ij}(x, y)$ . From (22), the $x$ component of the total relative atomic lattice displacement between the atomic planes at $y = \pm a/2$ reads:

$$\delta (x)\equiv u_z(x,+a\mathrm{/}2)-u_z(x,-a\mathrm{/}2)=\frac{a}{\mu }[{\bar{\sigma }}_a(x)+\frac{\mu b}{\pi a}\arctan \left(\frac{a}{2x}\right)]+b\theta (-x)\mathrm{.}(23)\backslash delta(x)\; \backslash equiv\; u\_z(x,\; +a/2)\; -\; u\_z(x,\; -a/2)\; =\; \backslash frac\{a\}\{\backslash mu\}\; \backslash left[\; \backslash bar\{\backslash sigma\}\_a(x)\; +\; \backslash frac\{\backslash mu\; b\}\{\backslash pi\; a\}\; \backslash arctan\; \backslash left(\; \backslash frac\{a\}\{2x\}\; \backslash right)\; \backslash right]\; +\; b\; \backslash theta(-x).\; \backslash quad\; (23)$$

Furthermore on account of (21) the average shear stress $\bar{\sigma}_{yz}(x)$ in the layer is:

$${\bar{\sigma }}_{yz}(x)={\bar{\sigma }}_a(x)+\frac{\mu b}{2\pi }\frac{1}{a}{\int }_{-a\mathrm{/}2}^{+a\mathrm{/}2}dy\frac{x}{x^2+y^2}={\bar{\sigma }}_a(x)+\frac{\mu b}{\pi a}\arctan \left(\frac{a}{2x}\right)\mathrm{.}(24)\backslash bar\{\backslash sigma\}\_\{yz\}(x)\; =\; \backslash bar\{\backslash sigma\}\_a(x)\; +\; \backslash frac\{\backslash mu\; b\}\{2\backslash pi\}\; \backslash frac\{1\}\{a\}\; \backslash int\_\{-a/2\}^\{+a/2\}\; dy\; \backslash frac\{x\}\{x^2\; +\; y^2\}\; =\; \backslash bar\{\backslash sigma\}\_a(x)\; +\; \backslash frac\{\backslash mu\; b\}\{\backslash pi\; a\}\; \backslash arctan\; \backslash left(\; \backslash frac\{a\}\{2x\}\; \backslash right).\; \backslash quad\; (24)$$Comparison of (24) and (23) shows that:

$$\delta (x)=(a\mathrm{/}\mu ){\bar{\sigma }}_{yz}(x)+b\theta (-x)\mathrm{.}(25)\backslash delta(x)\; =\; (a/\backslash mu)\backslash bar\{\backslash sigma\}\_\{yz\}(x)\; +\; b\backslash theta(-x).\; \backslash quad\; (25)$$

Consider next an Eshelby screw dislocation with an extended core described by a continuous function $\eta(x)$ . The distortion $\beta^*$ becomes: $\beta_{yz}^*(x, y) = \delta_D(y)\eta(x) = -\delta_D(y)\int_x^{+\infty} dx \eta'(x)$ ; the Volterra dislocation corresponds to the limiting case $\eta(x) = b\theta(-x)$ . Displacements, strains and stresses are obtained by convolution using $d\beta_{yz}^*(x, y) \equiv -\delta_D(y)\eta'(x)dx$ [10] as elementary distortions. The analogs of Eqs. (24,25) are:

$${\bar{\sigma }}_{yz}(x)={\bar{\sigma }}_a(x)-\frac{\mu }{\pi a}{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}d{x}^{\mathrm{\prime }}\arctan \left(\frac{a}{2(x-x^{\mathrm{\prime }})}\right){\eta }^{\mathrm{\prime }}(x^{\mathrm{\prime }})(26)\backslash bar\{\backslash sigma\}\_\{yz\}(x)\; =\; \backslash bar\{\backslash sigma\}\_a(x)\; -\; \backslash frac\{\backslash mu\}\{\backslash pi\; a\}\; \backslash int\_\{-\backslash infty\}^\{+\backslash infty\}\; dx\text{'}\; \backslash arctan\backslash left(\backslash frac\{a\}\{2(x-x\text{'})\}\backslash right)\; \backslash eta\text{'}(x\text{'})\; \backslash quad\; (26)$$

$$\delta (x)=(a\mathrm{/}\mu ){\bar{\sigma }}_{yz}(x)+\eta (x)\mathrm{.}(27)\backslash delta(x)\; =\; (a/\backslash mu)\backslash bar\{\backslash sigma\}\_\{yz\}(x)\; +\; \backslash eta(x).\; \backslash quad\; (27)$$

Equation (27), which we use in the paper, shows the relation between the coarse-grained displacement $\delta$ , and the discontinuity $\eta$ .

References

1. 1. Carpio, A. and Bonilla, L.L.: Discrete models for dislocations and their motion in cubic crystals. Phys. Rev. B 12, 2005, 1087–1097.
2. 2. Carstensen, C., Hackl K. and Mielke, A.: Nonconvex potentials and microstructures in finite-strain plasticity. Proc. R. Soc. Lond. A 458, 2002, 299–317.
3. 3. Choksi, R., Del Piero, G., Fonseca, I. and Owen, D.R.: Structural deformations as energy minimizers in models of fracture and hysteresis. Math. Mech. Solids 4, 1999, 321–356.
4. 4. Conti, S. and Zanzotto, G.: A variational model for reconstructive phase transformations, and their relation to dislocations and plasticity. Arch. Rational Mech. Anal. 173, 2004, 69–88.
5. 5. Christian, J.W. and Vitek, V.: Dislocations and stacking faults. Rep. Prog. Phys. 33, 1970, 307–411.
6. 6. Del Piero, G. and Truskinovsky L.: Macro and micro-cracking in 1D elasticity. Int. J. Solids Struct. 38, 2001, 1135–1148.
7. 7. Denoual, C.: Dynamic dislocation modeling by combining Peierls–Nabarro and Galerkin methods. Phys. Rev. B 70, 2004, 024106.
8. 8. Denoual, C.: Modeling dislocations by coupling Peierls–Nabarro and element-free Galerkin methods. Comput. Meth. Appl. Mech. Engrg. 196, 2007, 1915–1923.
9. 9. Ericksen, J.: Equilibrium of bars. J. Elast. 5, 1975, 191–202.
10. 10. Eshelby, J.D.: Uniformly moving dislocations. Proc. Phys. Soc. London A 62, 1949, 307–314.
11. 11. Hakim, V. and Karma, A.: Crack path prediction in anisotropic brittle materials. Phys. Rev. Lett. 95, 2005, 235501.
12. 12. Hirth, J.P. and Lothe J.: Theory of dislocations, 2nd edn. Wiley & Sons, New York, 1982.
13. 13. Miller, R., Phillips, R., Beltz, G. and Ortiz, M.: A non-local formulation of the Peierls dislocation model. J. Mech. Phys. Solids 46, 1998, 1845–1867.
14. 14. Mura, T.: Micromechanics of defects in solids, 2nd edn. Martinus Nijhof, Dordrecht, 1987.
15. 15. Nabarro, F.R.N.: Dislocations in a simple cubic lattice. Proc. Phys. Soc. 59, 1947, 256–272.
16. 16. Ortiz, M. and Phillips, R.: Nanomechanics of defects in solids. Adv. Appl. Mech. 36, 1999, 1–79.
17. 17. Pellegrini, Y.-P., Denoual C. and Truskinovsky, L., in preparation.1. 18. Ponte Castañeda, P. and Suquet, P.: Nonlinear composites. Adv. Appl. Mech. 34, 2002, 171–302.
18. 19. Rice, J.R.: Dislocation nucleation from a crack tip: an analysis based on the Peierls concept. J. Mech. Phys. Solids 40, 1992, 239–271.
19. 20. Rockafellar, T.: Convex analysis. Princeton University Press, Princeton, 1997.
20. 21. Rosakis, P.: Supersonic dislocation from an augmented Peierls model. Phys. Rev. Lett. 86, 2001, 95–98.
21. 22. Stillinger, F.H. and Weber, T.A.: Packing structures and transitions in liquids and solids. Science 225, 1984, 983–989.
22. 23. Sun, Y., Beltz, G.E. and Rice, J.R.: Estimates from atomic models of tension-shear coupling in dislocation nucleation from a crack tip. Mat. Sci. Eng. A 170, 1993, 69–85.
23. 24. Truskinovsky, L.: Kinks versus shocks. In: Fosdick, R., Dunn, E. and Slemrod, M. (Eds.), Shock Induced Transitions and Phase Structures in General Media, IMA, Vol. 52, Springer-Verlag, 1993, pp. 185–229.
24. 25. Truskinovsky, L.: Fracture as a phase transformation. In: Batra, R. and Beatty, M. (Eds.), Contemporary research in mechanics and mathematics of materials, CIMNE, Barcelona, 1996, pp. 322–332.
25. 26. Wang, Y. and Khachaturyan, A.G.: Three-dimensional field model and computer modeling of martensitic transformations. Acta Mater. 45, 1997, 759–773.
26. 27. Woodward, C.: First-principles simulations of dislocation cores. Mat. Sci. Engrg. A 400–401, 2005, 59–67.