Search tips
Search criteria 


Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
Proc Math Phys Eng Sci. 2016 April; 472(2188): 20160039.
PMCID: PMC4892284

Lattice continuum and diffusional creep


Diffusional creep is characterized by growth/disappearance of lattice planes at the crystal boundaries that serve as sources/sinks of vacancies, and by diffusion of vacancies. The lattice continuum theory developed here represents a natural and intuitive framework for the analysis of diffusion in crystals and lattice growth/loss at the boundaries. The formulation includes the definition of the Lagrangian reference configuration for the newly created lattice, the transport theorem and the definition of the creep rate tensor for a polycrystal as a piecewise uniform, discontinuous field. The values associated with each crystalline grain are related to the normal diffusional flux at grain boundaries. The governing equations for Nabarro–Herring creep are derived with coupled diffusion and elasticity with compositional eigenstrain. Both, bulk diffusional dissipation and boundary dissipation accompanying vacancy nucleation and absorption, are considered, but the latter is found to be negligible. For periodic arrangements of grains, diffusion formally decouples from elasticity but at the cost of a complicated boundary condition. The equilibrium of deviatorically stressed polycrystals is impossible without inclusion of interface energies. The secondary creep rate estimates correspond to the standard Nabarro–Herring model, and the volumetric creep is small. The initial (primary) creep rate is estimated to be much larger than the secondary creep rate.

Keywords: lattice growth, moving boundaries, vacancy diffusion, vacancy source, continuum kinematics

1. Introduction

At high homologous temperatures, polycrystalline solids exhibit time-dependent inelastic deformation—creep, even when subjected to stresses much lower than the nominal yield strength. Continuing creep eventually results in cavitation and rupture (figure 1a). The creep properties of materials are the major limiting factor in high temperature applications.

Figure 1.
(a) Macroscopic observations of creep in a uniaxial step loading experiment. The instantaneous elastic deformation is followed by the transient primary creep, then by the secondary creep with constant creep rate and finally the tertiary creep with rapidly ...

Depending on the temperature and stress level, several micromechanisms of creep can be distinguished [1,2] including

  •  Nabarro–Herring creep: diffusion of vacancies through the crystalline lattice and complementary diffusion of atoms (figure 1b) through the vacancy–atom exchange mechanism, leads to the creation of new lattice (figure 1c) layers on some boundaries and disappearance of lattice layers on other boundaries. Boundaries are sources and sinks of vacancies.
  •  Coble creep: diffusion of vacancies/atoms along grain boundaries with the same outcome: lattice growth/disappearance.
  •  Dislocation creep: vacancy diffusion enables dislocation climb and glide. A variety of mechanisms have been proposed, including the low stress Harper–Dorn creep.

Identification of the creep mechanisms, based on the macroscopic measurements, has sometimes proved difficult. The creep rate for dislocation-mediated Harper–Dorn creep exhibits the same linear dependence on stress as the diffusional creep rate, so that the only distinguishing feature between the Harper–Dorn, Coble and Nabarro–Herring creep mechanisms is the creep rate dependence on grain size [3]. In fact, the creep mechanisms are always to some extent speculative. On the one hand, experimental observations at high temperatures are macroscopic, and microscopic analysis is only available at the end of the experiment, after cooling. On the other hand, atomistic simulations of diffusion processes (rare events) in significant volumes are still mostly beyond the current state of the art. This underlines the importance of having a physically unambiguous and mathematically rigorous mesoscale continuum model for diffusional creep, capable of resolving the diffusion within the grain and the lattice growth/loss at the boundary.

It has long been recognized that a continuum theory of diffusion in crystals must take into account the crystalline lattice, and some elements of lattice kinematics have been used in the creep and diffusion models. Larché & Cahn [46] considered the interaction of elasticity and diffusion in great detail, but not the diffusional growth of the lattice. Berdichevsky et al. [7] recognized the contradiction between the kinematics of the usual mass continuum and the kinematics of diffusion through the lattice and attempted to resolve it by introducing ‘elastic velocity’. The resulting mathematical problem has similarities with the one developed here, but without elasticity–diffusion coupling. Garikipati et al. [8] attempted to define ‘lattice thermodynamics’. Without the proper lattice transport theorem, they failed to resolve the dichotomy between boundary kinematics and lattice kinematics, and were forced to adopt a diffuse interface model, with an interface energy field spread over some interface width; an approach followed later by Villani et al. [9].

With reference to figure 1b,c, the problem belongs to the class of moving boundary problems. A convenient formulation for numerical analysis of such problems is the phase field (diffuse interface) method. Such formulation has been proposed by Mishin et al. [10]. However, the mesoscale phase field models contain some arbitrariness unless they are firmly grounded in the mathematical structure of the sharp interface formulation, i.e. unless they converge to the well-defined sharp interface limit (in the sense of Gibbs excess quantities [11]) as the interface thickness vanishes. The full sharp interface formulation has remained elusive, partly because of uncertainty about the physical processes that governs the vacancy production/absorption at the boundaries, but also because of the incompleteness of the lattice continuum formulation.

In this paper, we first develop the full lattice continuum kinematics. Then, we formulate the problem of isothermal Nabarro–Herring creep with the goal of understanding the mathematical structure of the problem. Although simple analytic solutions are not expected, a complete formulation and dimensional analysis will provide some novel insights. The paper is organized as follows.

In the remainder of this section, we review the basic elements of standard continuum kinematics. In §2, we formulate the kinematics of lattice and interfaces, including

  •  the lattice (material) derivative,
  •  the transport theorem and mass balance,
  •  the definition of Lagrangian reference configuration for the newly created lattice, and
  •  the proper definition of the creep deformation rate, qualitatively different from other components of deformation.

Thus defined lattice kinematics enables the formulation of power balance for isothermal Nabarro–Herring creep. In §3, two dissipative mechanisms are identified: the bulk diffusion and the vacancy creation at the boundary. The governing equations are derived in §4 by means of the principle of virtual power. The coupled elasticity–diffusion problem is formally decoupled in §6 for the case of periodic grain assemblies, to clarify the mathematical structure of the standalone diffusion problem. In §7, dimensional analysis and simple approximate scaling are used to estimate the creep rate, and to further simplify the governing equations for the problem without boundary dissipation. The role of interface energy is analysed in §8, followed by summary and conclusions in §9.

(a) Review of the mass continuum

For reasons that will be apparent shortly, the standard continuum models can be called mass continuum. The quantity associated with the material is the mass density, i.e. this is the quantity advected by the material velocity v¯. Hence, for any material field Y (x,t), where x represents the spatial (Eulerian) coordinates, the material derivative (the time derivative following the material, i.e. the mass density) is defined as


Some variables and operators are different in mass and lattice continuum formulations. This difference is emphasized by using the overbar for the mass continuum formulation. The local mass conservation is expressed as


The special role of the mass density among all material fields is emphasized by the simple form of the transport theorem. For a material field Y (x,t) (i.e. given per unit mass)


where V¯(t) is the material volume, and its dependence on time signifies that the material occupies different spatial volume at different times.

In the fluid mix of multiple components, the natural choice for the material velocity is the barycentric velocity [12], which is defined from the requirement that the total momentum density ρv¯ is equal to the sum of the components’ momentum densities,1 thus ensuring that the linear momentum balance has the well-known Cauchy form (without body forces)


where σ¯ is the Cauchy stress. The simple expression for the rate of change of linear momentum, the right-hand side of (1.4), is the main reason for identification of material with mass density.

Diffusional creep is typically a very slow process and the momentum and its rate of change can be neglected. In such quasi-static (quasi-equilibrium) processes, kinetic energy and inertial forces are negligible. In other words, the inertial mass becomes irrelevant and only the mass as quantity of matter remains. Thus, the variable representing the material need not be the mass density. On the contrary, for the mathematical description of diffusion through the lattice, the choice of mass density as the quantity advected as material is distinctly inconvenient and counterintuitive.

In this paper, we replace the mass density with the lattice site density as the main material field, and develop the lattice continuum formulation for quasi-static processes, with all the essential continuum theorems. We emphasize that downgrading the role of mass density does not imply that mass conservation and transport are ignored; only the inertial aspect of mass is neglected.

2. Lattice kinematics of a crystal and its boundaries

We begin the analysis with a finite single crystal with a conserved mass, and we consider isothermal processes.2 Consider a single species of atoms and vacancies, without interstitials.3 At any instant t, the lattice velocity field is defined as suitable interpolation between the velocities of the lattice points. It is an Eulerian variable v(x,t), defined in the open domain Ω, its closure [partial differential]Ω having the outer normal n4. The lattice velocity gradient L is defined as


The velocity field characterizing the continuum carries with itself the lattice density NL, so that, in contrast to the standard mass continuum and its material derivative, we define the lattice derivative as the time derivative following the lattice. For an arbitrary field Y (x,t), which follows the lattice


The local conservation of lattice density is expressed as


In addition, the lattice density NL, the lattice point5 currently located at x, carries the vacancy concentration c (the fraction of lattice sites occupied by vacancies), and the deformation gradient F, which evolves as


At this point, it is instructive to consider the physical interpretation of deformation and velocity gradients, and provide a more rigorous definition of deformation of a material (lattice) element. We think of lattices as a networks (graphs) consisting of lattice points and the lines connecting those points. The lines may be based on nearest neighbour tessellation (Delaunay graph) or interatomic bonding. In general, the lattice graph is not unique; the interatomic bonds may be picked selectively, and Delaunay tessellation of a periodic set of points is not unique. Nevertheless, once defined, the graph has a distinct topology based on the connectivity of the network. This topology is invariant under the change of species that occupies any particular lattice site.

(a) Plasticity and phase transformations

The deformation gradient can be decomposed into the elastic-compositional deformation pseudo-gradient, Fec, and the plastic deformation pseudo-gradient, Fp. The former represents the motion which produces a deformed lattice topologically equivalent to the reference lattice, the latter changes the topology, i.e. the connectivity. The decomposition is multiplicative and analogous to the standard elastic–plastic decomposition [13]


The two components are in general incompatible (i.e. neither is a true gradient of a single-valued vector field), but the total gradient is compatible. The corresponding decomposition of lattice velocity gradient follows the standard crystal plasticity6 [14]


The physical content of the two deformation pseudo-gradients can be determined based on the above definition. The elastic-compositional deformation pseudo-gradient, Fec, includes lattice stretching and distortion. Such stretching and distortion of an infinitesimal volume element may be the result of

  •  Elastic deformation resulting from the action of stresses without change in composition (lattice site occupancy).
  •  Compositional deformationof a stress-free element resulting from changes in composition; in the present case, the vacancy concentration.

The plastic deformation pseudo-gradient, Fp, incorporates the changes in lattice topology, including the changes caused by dislocation glide and climb. Phase transformations can belong to either of the two pseudo-gradients, depending on whether they are accompanied by changes in topology (martensitic) or not.

We emphasize that the above deformation analysis is local; it only describes deformation at a material point within a crystal. If an assembly of grains (a polycrystal) is considered, the global deformation of the assembly caused by growth/disappearance of lattice sites at various parts of the grain boundary network, must be included. This portion of the total deformation cannot be described locally. We will use the word creep for this portion of deformation.

(b) Reference configuration

To keep the ensuing thermodynamic analysis simple, we consider only the elastic-compositional deformation gradient, i.e. we do not consider dislocation motion and phase transformations


To extend the definition of lattice-advected variables to the portion of the lattice currently growing at the boundary, it is sufficient to assert the continuity of the deformation gradient and the vacancy concentration within the crystal. Then, the growing portion of the lattice is characterized by the current boundary values of deformation gradient and vacancy concentration at the point of growth. Deformation gradient in a newly grown lattice being thus defined, the reference position X of a newly grown lattice point is formally defined by integrating7


Let NL0 be the number of lattice point per unit volume of the undeformed lattice.8 Then, the number of lattice point per unit volume of the deformed lattice is


The number of atoms per unit deformed volume is (1−c)NL.

(c) Boundary kinematics and creep deformation

The boundaries of the crystal serve as vacancy sources and sinks. The vacancy diffusion takes place by the vacancy–atom exchange mechanism and is therefore accompanied by the complementary diffusion of atoms. The boundaries cannot absorb or produce atoms. Instead, the lattice grows at the parts of the boundary where atoms arrive from the bulk of the crystal (vacancy source) and vanishes where atoms leave towards the bulk of the crystal (vacancy sink).

We consider the boundary velocity V as a field defined on the closure set [partial differential]Ω, and independent of the lattice velocity field. The difference in normal components of the two velocities is the lattice growth rate


Tentatively,9 we define the tangential component of boundary velocity to coincide with the lattice velocity component


where nn is a dyadic product.

The creep velocity gradient is the rank-2 tensor, constant within the domain Ω and defined on the basis of the difference between the boundary velocity V, and the lattice velocity at the boundary10


where left angle bracketLright angle bracketΩ is the average lattice velocity gradient in the crystal domain Ω. The creep deformation rate tensor is the symmetric portion of (2.12). For the simple definition of boundary velocity (2.11), the creep deformation rate is


Note that in contrast to the locally defined elastic-compositional deformation rate, the creep deformation rate is piecewise constant and discontinuous across grain boundaries. This is not unexpected and it conforms to the physical understanding of diffusional creep.

(d) Diffusion

In the mix of vacancies and atoms with mass m, the mass density is given as


so that, following (2.2),


Because vacancies are massless, the average velocity of atoms va is equal to the barycentric velocity of the mix. The local mass conservation requires


Upon comparing the last two results


The vacancy flux J, relative to the lattice, is the negative of the atomic flux Ja. It is a vector advected with the lattice, and can be equivalently defined from relative velocities of either atoms or vacancies


where vv is the average velocity of vacancies. Then, we obtain the diffusion equation as


(e) Transport theorem and mass conservation

Let the quantity Y be given per lattice site. Consider the lattice volume V (t), i.e. the volume that follows the prescribed set of lattice sites. The lattice continuum variant of the transport theorem is


For a variable crystal domain Ω(t), with lattice growing/vanishing at the boundary [partial differential]Ω(t), the transport theorem takes the form


Using (2.11), (2.14), (2.19), the mass of a crystal grain will change as


If a single grain is to satisfy mass conservation,11 then, upon substitution of (2.9) and (2.17) and application of the divergence theorem, we obtain


The integral condition (2.21) does not prevent local atom exchange between the grains, nor the grain boundary diffusion, as long as the net effect of those on the grain mass vanishes. Without a model for grain boundary diffusion, the condition (2.21) allows the grain boundary to be a vacancy reservoir, with instantaneous redistribution of vacancies around the grain. The stronger, local condition,


implies absence of boundary diffusion by requiring that lattice growth/disappearance at the boundary be directly related to the normal atomic flux. In fact, the global condition (2.21) is insufficient to define the problem; additional information about redistribution of vacancies on the boundary is required. For now, we will consider the simpler, local condition (2.22). The creep deformation rate can then be written as


3. Energy (rate) balance for a single crystal

(a) Free energy

At constant temperature, the energy density can be represented as function of concentration and deformation gradient, Φ(F,c), given per unit reference volume, 0=dΩNL/NL0. The total free energy12 of the grain with the current volume Ω is


Using (2.19), the rate of change of the free energy of a grain is


The first term can be expressed as


where [partial differential]Φ/[partial differential]F is the lattice nominal stress, and σ is the lattice Cauchy stress13


With the shorthand: Φ′=[partial differential]Φ/[partial differential]c, the rate of change of the grain free energy is


(b) Dissipation and the second law

Local diffusional dissipation rate in the lattice is modelled simply as proportional to the local rate of change of concentration, so that the total bulk dissipation rate is


For convenience, the diffusion potential M is defined with respect to the reference configuration. In addition to the bulk dissipation, the creation and absorption of vacancies at the crystal boundaries must be a dissipative process. The physical mechanisms of boundary dissipation will be discussed later in this section. For now, we note that, following (2.22), the dissipation rate at a boundary point can be equivalently assumed to be proportional either to the lattice growth rate g, or to the normal vacancy flux Jn


The total dissipation rate is then


For isothermal processes, the second law of thermodynamics requires


The linear constitutive law with positive mobility B,


guarantees that the volume portion of (3.8) is non-negative. Similarly, the boundary integral in (3.8) will be non-negative if


In SI units, the dimensions of mobilities B and b are m4(Ns)−1 and m3(Ns)−1, respectively.

The boundary dissipation has the usual Rayleigh dissipation structure with the potential D(Jnv), and the complementary (numerically equal) potential D¯(χ+M)


so that


(c) Physical mechanism of boundary dissipation

In figure 2, two examples of relaxed grain boundaries are shown. The relaxed grain boundaries can be characterized by the excess volume, in the sense of Gibbs surface excess quantities [11]. It is the excess volume per atom in the interface layers with respect to the volume per atom in the bulk. The surface excess volume ΔV is the conjugate of pressure with respect to the surface energy γ [17]:


where T is the absolute temperature. Moreover, it correlates with the interface energy [18]. When comparing the grain boundaries in a single material, larger excess volume implies larger interface energy.

Figure 2.
Examples of relaxed asymmetric tilt Σ9left angle bracket110right angle bracket grain boundaries in Cu, with moderate (a) and large (b) excess volume. Colour scale is based on centrosymmetry parameter [16]. (Online version in colour.)

Imagine now that interfaces in figure 2 are stretched under tensile stresses, but in such a way that topology is preserved (i.e. the nearest neighbours remain the same). The excess volume increases, so that the local energy of atoms is high and can be lowered by a lattice atom popping into the large open space at the boundary, such as those indicated by arrows in figure 2b. Because the moving atom leaves a vacancy behind, this amounts to vacancy nucleation. The event is thermally activated and thus dissipative.

Consider the thermally activated process, with the activation energy for the creation/absorption of a vacancy Qv, the temperature T, the Boltzmann constant k and the attempt frequency n0. Let the lattice parameter be N−1/3L=λ. Assuming that all newly created vacancies are transported towards the bulk, the normal vacancy flux is proportional to the rate of creation of vacancies, n. Assuming further that the density of potential nucleation sites at the interface is approximately one, we obtain


Similarly, assuming that the concentration is close to the equilibrium value for the given temperature, c0, the bulk mobility can be written as


where we assume the same attempt frequency n0, Qd is the activation energy for elementary atom–vacancy exchange event, and the density of potential atom–vacancy exchange sites is c0. For later use, we note that


We emphasize that the thermal activation process considered here is limited to the nucleation of a vacancy which is immediately transported into the crystal. Nucleation of a stable vacancy that remains at the boundary would have much larger activation energy.

(d) Power balance

Let the crystal be loaded on the boundary with tractions t. The power balance for an isothermal process can be written as


From (2.10), (2.11) and (2.22), we obtain


Then, with (3.5, 3.6), the power balance (3.17) takes form


4. Virtual power and the governing equations for a single crystal

At any instant t, the fields v(x,t) and J(x,t) in (3.19) are subject to independent variations δv(x,t) and δJ(x,t). The principle of virtual power is written as


The usual manipulation and the argument based on independence of variations yield four independent variational statements


which, in turn, yield the governing equations. The first pair yields the stress (quasi-) equilibrium with boundary conditions:


The second pair implies


With thus defined diffusion potential M and the boundary flux conjugate χ, it is instructive to re-examine the requirement that the boundary dissipation is non-negative,


We recall that the work of boundary tractions (3.13) can be represented as the sum of the work on lattice deformation, t[center dot]v dt, and the work on lattice creation, tng dt. With (4.4) and (2.22), the boundary dissipation condition (4.5) can then be re-written as


The difference between the left- and right-hand sides of this inequality is the boundary dissipation associated with the vacancy nucleation/absorption.

The diffusion is then governed by


The boundary condition in (4.7) also determines the lattice growth rate (2.22). Upon expressing the normal flux from (3.10), we obtain


5. Partially linearized problem

While the variable domain represents an irreducible nonlinearity, the problem can be simplified for small elastic-compositional strains and small vacancy concentration. Moreover, we choose the coordinate system at the initial centre of mass of the grain and ignore the rigid body motion, implying small rotations. The elastic-compositional deformation gradient can be expressed as


For small elastic deformation, we take Φ(ε,c). Following Larché & Cahn [46], we express the total strain ε as the sum of elastic strain and compositional volumetric strain: −η(cc0), (η>0 for vacancies), the free energy density may be expressed as the sum of elastic energy and mixing energy


In (5.2), c0 is the equilibrium concentration (for unstressed crystal at a given temperature), the elastic strain is the difference between the total strain and the compositional strain: (ε+η(cc0)I), and E is the elastic stiffness tensor. Thus, the compositional strain plays the role of the eigenstrain (stress-free strain analogous to thermal dilation/contraction). Further,


The free energy of mixing, ψ(cc0), ψ(0)=0, can be expressed as the sum of bonding (atom–atom versus atom–vacancy) and entropic parts, as in the regular solution model [19]. We expand the energy of mixing about c=c0 up to the second order (appendix A), to obtain


With the Boltzmann constant k=1.38×10−23 Nm K−1, the temperature, T=103 K, typical lattice density for metals, NL0=8×1028 m−3 and the equilibrium vacancy concentration c0=10−3, we obtain κ≈1 TPa. The coefficient κ can be understood as the energetic penalty for departing from the equilibrium concentration at a given temperature. Based on its magnitude, we expect very small departures from the equilibrium concentration, as we assumed in connection to (3.15).

We note that ϕ/ϕ′ is of the order of elastic-compositional strain, whereas ψ/ψ′ is of the order (cc0). Assuming a small lattice velocity (v[center dot][nabla]c[double less-than sign][partial differential]c/[partial differential]t; n[center dot]v[double less-than sign]V n), the diffusion governing equations take the form


Finally, to complete the mathematical formulation of the problem, we express the total local strain (elastic+compositional) in terms of the lattice displacement


Here, we recall that the strain ε does not include the creep strain. Thus, we can assume that this strain is small. We can ensure that the displacements u(X) are small by subtracting the rigid body motion of the crystal. We emphasize that the displacements u(X) are defined even for the part of the lattice which was created during the creep deformation, by means of the inverse mapping (2.8).

Then, at each instant t, we have a linear elasticity problem with concentration dependent eigenstrain and external loading:


The diffusion–elasticity problem with moving domain boundaries exhibits coupling in both, differential equations and boundary conditions. In §6, we consider the special case of periodic polycrystals, whereby formal decoupling of diffusion and elasticity is possible, but at a cost of introducing non-locality into the diffusion boundary condition.

6. Periodic polycrystals

For periodic assemblies of grains, the instantaneous linear elasticity eigenstrain problem can be solved using the methods developed for periodic microstructures such that the geometric centres of the grains form a Bravais lattice [2022]. First, decompose the elasticity problem into two problems with different characteristic length scales with microscale displacement uc and macroscale displacement u0. If the uniform stresses σ0 are prescribed, the elasticity problem decouples from the diffusion problem. The microscale elasticity problem is the eigenstrain problem whereby the averages of displacement uc and the corresponding stress σc, over the domain Ω, vanish.

(a) Solution to the elasticity problem

We consider elastically isotropic grains:


where λ and μ are the Lamé elastic constants. Let


and let the symbol left angle bracket[center dot]right angle bracket denote the lattice volume average. For example,


Under uniform composition, the stress-free assembly undergoes uniform compositional volumetric strain


Define the relative concentration as


For a given periodic field c~(x), with vanishing volume average over the domain Ω, the displacement uc is also periodic, and the compositional stresses are given as


The asterisk symbol represents the convolution


The kernel Γij(xξ) is given in appendix B in the form of Fourier series on the Bravais lattice. It is evident from (5.7) that the compositional stresses must be homogeneous linear functionals of c¯ over the domain Ω; the uniform c¯=c¯ with no applied tractions produces no stresses. The compositional pressure (derived in appendix B) exhibits a local dependence on c~


The volume average of the compositional pressure vanishes.

(b) Decoupled diffusion problem

For the problem with applied far-field traction, the stress σ0 is uniform. We can eliminate the compositional pressure from the governing equations (5.5)


We recall that the compositional traction tnc=nσc hides the volume convolution (6.7), whereas c¯ is the volume average. Thus, the problem is described by the standard diffusion differential equation, but on a domain with moving boundaries and with integrodifferential boundary conditions. If the grains remain polyhedral and faces do not rotate, then the normal velocity V n is uniform at each face, so that that the normal gradient of concentration must be uniform at each face.

(c) Equilibrium

Consider now the equilibrium under uniform applied stress. The concentration must be uniform: c¯eq=c¯eq. Denote the tractions arising from deviatoric stresses with a prime


The equilibrium implies


This condition must be satisfied on each (flat) face of a crystal. This is clearly possible only if the applied stress is hydrostatic σ0=−p0I in that case, the equilibrium concentration is given as


As expected, the equilibrium vacancy concentration is pressure-dependent.

Under a deviatoric loading, σ0≠−p0I, no equilibrium is possible. Thus, a deviatorically stressed crystal can never be in thermodynamic equilibrium with its boundaries14 (unless interface energies are introduced, which we consider in §8).

The above conclusion is in apparent contradiction to the analysis of Larché & Cahn [5], who concluded that such equilibrium is possible. The source of this apparent paradox is as follows. While Larché & Cahn [5] acknowledge the possibility of interface motion relative to the lattice, they only consider lattice growth/loss as an exchange between two neighbouring crystals (see §3.4 of their paper). Most importantly, their variational statement (equation (11) of the 1978 paper) does not include lattice growth/loss. Specifically, the term analogous to the fourth equation of (4.2) is missing in their derivation.15

(d) Volumetric creep

Consider the uniaxial step loading of a polycrystal (figure 1a). As vacancies are nucleated at boundaries under tension, they propagate into the crystal. The equilibrium concentration for the given pressure (6.12) indicates that some creep volume change is inevitable. From (6.4) and (6.12), and for uniaxial applied stress σ0, the expected total creep volume change is


which, with the earlier estimate for κ (5.4) and η of the order 10−1, is smaller than or, at most, equal in magnitude to the elastic strains.

7. Dimensional analysis and scaling

The governing equation (6.9) represents a problem with moving boundaries and complex integrodifferential boundary conditions and thus hold little hope for an analytic solution even for the simple rectilinear geometry. A systematic numerical treatment, beyond the scope of this paper, is needed. Nevertheless, dimensional analysis can provide some insights.

(a) Characteristic times

Let the characteristic grain size be a0. The characteristic diffusion time scale is


With t^=t/τdiff and ^=a0, the non-dimensional equations take the form


Another characteristic time, associated with the vacancy nucleation/absorption on boundaries, emerges from the analysis


With (3.16), we find


(b) Simplified equations in the absence of boundary dissipation

We expect that for close-packed metals QdQv, so that the τdiff/τvac is very large. For covalently bonded ceramics or with impurities segregated on the boundaries, it is possible that Qv>Qd. But even then, the ratio τdiff/τvac can be comparable to unity only for very small crystals and very large values of Qv. Thus, for most cases of interest, the boundary produces/absorbs as many vacancies as the lattice can transport, so that the expression in brackets in (7.2) must be very small. Requiring that it vanishes amounts to explicitly neglecting the boundary dissipation, so that the simplified diffusion problem16 is


with the boundary moving as


We note that the compositional tractions scale as tncηβc¯, so that their effect is typically not negligible. Moreover, it is evident from (7.5) that a non-vanishing solution with continuous boundary concentration is possible only with the compositional tractions. Without them, the concentration must have discontinuities at the edges between the faces. Such solutions are possible for the diffusion equation [23].

(c) Secondary creep rate

To estimate the rate of creep, we consider the simple periodic geometry with grains being orthogonal parallelepipeds, as shown in figure 3. The assembly is loaded with the applied stress field σ10, and the characteristic grain size is a0=(a1a22)1/3. We note that, in this configuration, the average of the compositional tractions tnc over each face Sk must vanish: tnck=0. Thus, in the boundary condition (7.5), the fluctuations of the boundary concentration around its face average c¯k, multiplied by κ~, cancel out the fluctuations in the compositional tractions. The face-average boundary condition then reads


so that, from (7.6), the secondary creep rate is approximately


Assume that in the secondary creep the volume average of concentration has reached its equilibrium value (6.12). Upon subtracting the boundary condition on S2 from the one on S1, we obtain


The concentration gradient is expected to scale as


so that the creep rate is


which is the usual estimate for Nabarro–Herring creep [1].

Figure 3.
Crystalline grain in the shape of orthogonal parallelepiped, loaded in uniaxial tension. Notation used in (7.7)–(7.11) is shown.

(d) Primary creep

Consider the configuration in figure 3, with initially cubic grains and with the uniaxial step loading at t=0. The instantaneous response is characterized by the instantaneous source solution [23], singular at t=0. Theoretically, the initial boundary velocity is infinite, but only for an instant, which is not measurable.17 A more productive analysis follows from consideration of the boundary condition (7.5). For a short period of time, the boundary concentration at face S1 (except close to the edges) must be


While most of the domain is still at uniform c0, the initial boundary flux may be estimated by taking (for example) the lattice parameter λ as the characteristic length


The initial creep rate is larger than the secondary creep rate and is proportional to 1/a0:


This author is not aware of systematic experimental measurements of the primary creep rates that would confirm or falsify (7.14), but such tests of grain size dependence are certainly feasible.

8. Effects of interface energy

(a) Equilibrium

In the present formulation, the interface energy is not included in the total energy balance (§3). Its effects on creep are expected to be significant only for small grains. However, one interesting effect independent of the grain size is that it makes equilibrium of a deviatorically stressed grain assemblies possible. As the aspect ratio of grains increases, the continuing creep produces new interface area at an increasing rate. At a certain point, the work provided by the applied stress becomes insufficient to compensate for the increasing interface area.

To formally derive the governing equations with interface energy γ, additional terms would be added in the power balance (3.19). To illustrate the general principle, we consider the assembly of orthogonal parallelepipeds, subjected to uniaxial stress σ10, as shown in figure 3. The power supplied by the applied stress balances the rate of surface energy increase when


The factor 1/2 splits the interface energy between the two adjacent grains. Assuming that the creep deformation is purely deviatoric, we derive the equilibrium condition


The relevant stress levels for diffusional flow are of the order 10−6 to 10−5 μ≈0.2−2 MPa [1]. Consider the aspect ratio a1/a2=2, and typical interface energies of the order of 1 N m−1. For the critical stress to be in the range relevant for diffusional creep, the grain size would have to be a2=2.5–0.25 μm. Thus, an initially cubic grain of size a0=1.26 μm with interface energy of 1 N m−1, subjected to uniaxial tension (figure 3) will elongate by creep until a1=2 μm; a2=1 μm; a1/a0=1.59. After the primary creep, the creep rate should be diminishing gradually as a result the interface energy effect. However, this is not observed in the typical creep experiment. Instead, the transition to tertiary creep is observed with accelerated creep, cavitation and eventually—rupture.

By computing the stresses inside the grain produced by surface tension and their inclusion in (6.11) or (7.5), it is easily verified that this equilibrium is in fact characterized by hydrostatic stress state, as the deviatoric stresses are cancelled out by the surface tension effects.18 If we envision the polycrystal as the closed network of boundary walls with tension in walls equilibrated by compression of crystals inside the cells, the total stress tensor in crystal (under far-field tractions) can be written as the sum of applied, compositional and surface tension stresses σγ. The equilibrium condition (uniform diffusion potential) can be written as the requirement that the tractions on each face are the same


Thus, the theoretical arrest of creep is the result of the diminishing diffusion capacity of the lattice as its stress state approaches the hydrostatic state. However, the interfaces which were originally under tension, such as S1 in figure 3, are still under undiminished tension, which promotes vacancy nucleation. This may be the contributing factor in cavitation and rupture which characterizes the tertiary creep. Instead of being absorbed by the lattice, nucleated vacancies congregate at the interface, causing stress concentration and continuing macroscopic creep. It should be noted that the present analysis of the vacancy nucleation at a boundary (§3) is based on the assumption of its immediate absorption by the lattice. The energetics of stationary vacancy nucleation at the boundary is very different. For such purpose, the dependence of interface energy on the vacancy content (excess volume) must be defined.

(b) Flat boundaries

While the effect of surface energy on the stress state in a grain is important only for small elongated grains, its effect in promoting flat interfaces between crystals is ubiquitous. The grain boundaries typically form with low energy orientation characterized by cusp-like minima [17], which makes curving of interfaces under conditions of diffusional creep energetically unfavourable.

9. Summary and discussion

The lattice continuum theory developed in this paper, represent natural and intuitive framework for analysis of diffusion in crystals and lattice growth/loss at the boundaries. The formulation includes the definition of the Lagrangian reference configuration for the newly created lattice, transport theorem and mass conservation. The latter connects the creep rate tensor, constant within a crystal, to the normal diffusional flux at the boundaries. The linearized version of the creep deformation rate (2.23) can be written as


The kinematics is extended to include plastic deformation of the crystal (§2), although this problem is not pursued in the remainder of the paper.

Within this framework, the governing equations for Nabarro–Herring creep under isothermal conditions are developed with coupled diffusion and elasticity. At first, we considered both, the diffusional dissipation in the bulk, and the boundary dissipation accompanying vacancy nucleation/absorption. Upon concluding that the boundary dissipation is usually negligible, the simplified diffusion problem is characterized by the boundary condition equating the diffusion potential (equal everywhere to the rate of change of energy density with concentration) to the normal boundary traction. The linearized version reads


Although formally decoupled from elasticity for the case of periodic arrangement of grains, the boundary conditions (9.2), interpreted as conditions on boundary vacancy concentration, represent a set of integral equations involving the values of concentration throughout the volume. The condition (9.2), can also be interpreted as the requirement that the power of tractions expanded on lattice growth, is expanded locally on diffusion at the boundary. In the linearized version


While the problem with moving boundaries did not yield analytical solutions, the analysis and decoupling of governing equations and subsequent dimensional analysis provide some insights, including the following

  •  The secondary creep rate estimate corresponds to the standard Nabarro–Herring model.
  •  The volumetric creep is small.
  •  Equilibrium under deviatoric stresses is not possible without the inclusion of surface energy and then, only for grains with high aspect ratio.
  •  The main features of the primary creep are qualitatively understood. Future numerical simulations should provide quantification.

The list of unresolved issues and directions for future research include

  •  Inclusion of surface energy and surface stresses.
  •  Inclusion of boundary diffusion for Coble creep.
  •  While the growth of cavities is well understood [25,26], the conditions for their nucleation (i.e. the transition from secondary to tertiary creep) are not. For the flat faces that remain flat and do not rotate, the flux at each face must remain uniform. For continuing secondary creep without volume change C22=−C11/2, the normal inward flux of the face S1 in figure 3 must continuously increase: (−J1=C11a1), and so does the peak concentration in the centre of S1, together with the maximum local compositional tension [−pc in (6.8)]. This is qualitatively confirmed by Villani et al. [9] simulations. With the increasing concentration, the probability of cavity nucleation increases, leading to tertiary creep, growth of cavities and rupture.
  •  Related to the above issues is the question of vacancy nucleation without subsequent absorption by the lattice. The activation energy of such process will be much higher than the one characterizing the nucleation with lattice absorption. A dependence of surface energy of the vacancy content (excess volume) appears to be the key component of such process.
  •  Formulation of dislocation creep. While the bulk kinematics of plastic deformation is easily accommodated within the lattice continuum (§2), some questions remain open. In the bulk: description of edge dislocations as vacancy sinks and their climb. At the boundaries: continuity of state variables at the growing boundary, and dislocation-boundary reactions.


S. Forest provided helpful discussions and suggestions.

Appendix A. The free energy of mixing

Following the regular solution model [19], the energy of mixing (required to vanish at stress-free equilibrium) can be written

A 1

where k is the Boltzmann constant, T is the temperature and α is the parameter which depends on the nearest neighbour interactions and the crystal structure. The second term in (A 1) is the entropy of mixing. Then,


Finally, we obtain

A 2

Appendix B. Elasticity solution for periodic eigenstrain with the vanishing average

Consider an infinite periodic assembly of elastically isotropic grains subject to periodic volumetric eigenstrain (−ηc), with the vanishing volume average. With (6.1), the governing equations (5.7) are

B 1

The grains form a Bravais lattice. Let the vectors of the reciprocal lattice be denoted K. The normalized concentration is a periodic function on the Bravais lattice and can be expanded into a three-dimensional Fourier series

B 2

The components of the composition gradient are

B 3

Displacement is also a periodic function on the same Bravais lattice

B 4

Upon substituting (B 2)–(B 4) into (B 3), we obtain for each K

B 5

where Kj are Cartesian components of K, and K=|K|. The system of equations (B 5) is readily solved

B 6

Then, we compute the stresses arising from the eigenstrain distribution as

B 7

or, with the definition of γ(K) (B 2),

B 8

with the convolution operator

B 9

We note that Γii(xξ)=δ(xξ), so that the compositional pressure is given in the local form as

B 10


1The kinetic energy density of the mix does not have the simple form ρv¯2/2, but this is usually ignored in the dynamic theories, with the justification that the kinetic energy of diffusion is negligible.

2Most (if not all) parameters are temperature dependent. Coupling of heat conduction with the present formulation does not appear to present a conceptual problem.

3To minimize the book-keeping, we consider only two specific substitutional species. The formulation is readily generalized to multiple substitutional and interstitial species in a manner of Larché & Cahn [46]. Only the question of segregation of interstitial at the boundary requires further discussion.

4The growth/disappearance of the lattice on the closure [partial differential]Ω requires that a different velocity field, the boundary velocity, be defined on the closure. This will be addressed later in this section. In general, it is useful to think of all fields associated with the lattice as being defined on an open domain.

5Considered as continuum of points; discrete lattice only serves as an aid to intuition.

6Because the dislocation motion is defined with respect to the lattice, the crystal plasticity continua are, in fact, the lattice continua. In the absence of diffusion lattice and mass continua are identical, so that this distinction is usually not emphasized. However, the decomposition of mass velocity gradient analogous to (2.6) would be very different.

7This question requires more discussion if plasticity is included, particularly if the description of such deformation includes additional state variables. Continuity of state variables in a newly created lattice may not always be physically justified.

8The undeformed lattice is stress-free and with the reference vacancy concentration. The reference vacancy concentration is most conveniently chosen as the stress-free equilibrium concentration for a given temperature.

9Boundary sliding will not be considered here.

10Cf. Mesarovic & Padbidri [15] for the motivation for (2.12).

11Which need not be the case if grain growth is relevant, as in the problem of dynamic recrystallization.

12Interface energy plays a marginal role in the Nabarro–Herring creep. This will be quantified in §8.

13Nominal and Cauchy stress are defined as power conjugates to the corresponding kinematic rates, and are therefore different in mass and lattice continua.

14The rate of diffusion/creep may of course be immeasurably slow, so that the lack of equilibrium may not have any practical consequences.

15We note that this term is the result of the proper formulation of the lattice continuum transport theorem (2.19) and its application (3.2).

16The governing equations (7.5)–(7.6) for the problem without boundary dissipation can be derived directly by applying the principle of virtual power, as outlined in §4, but with the mass conservation (2.22) enforced by means of a Lagrange multiplier. In this case, the relative boundary motion is formally considered as an independent variable with variation δg, in addition to the flux and lattice velocity fields.

17This singularity should be interpreted with caution. It is in the form of instantaneous Dirac delta function on the boundary and it disappears immediately. In any case, the concentration is only defined in the open domain Ω; its boundary value is the limit as the boundary is approached from the inside of the crystal.

18In this qualitative discussion, we ignore the differences between surface stresses and surface energy in solids (cf. [24]).

Data accessibility

There are no additional data.

Competing interests

I have no competing interests.


The work was supported in parts by: visiting professorship at Centre des Matériaux, Mines ParisTech, France; US scholar grant by Fulbright Foundation; and, the professional leave programme at Washington State University.


1. Frost HJ, Asby MF 1982. Deformation-mechanism maps: the plasticity and creep of metals and ceramics. Oxford, UK: Pergamon.
2. Stouffer DC, Dame LT 1996. Inelastic deformation of metals. Hoboken, NJ: John Wiley & Sons.
3. Langdon TG. 2000. Indentifying creep mechanisms at low stresses. Mater. Sci. Eng. A 283, 266–273. (doi:10.1016/S0921-5093(00)00624-9)
4. Larché FC, Cahn JW 1973. A linear theory of thermochemical equilibrium of solids under stress. Acta Metall. 21, 1051–1063. (doi:10.1016/0001-6160(73)90021-7)
5. Larché FC, Cahn JW 1978. Thermochemical equilibrium of multiphase solids under stress. Acta Metall. 26, 1579–1589. (doi:10.1016/0001-6160(78)90067-6)
6. Larché FC, Cahn JW 1985. Interaction of composition and stress in crystalline solids. Acta Metall. 33, 331–357. (doi:10.1016/0001-6160(85)90077-X)
7. Berdichevsky V, Hazzledine P, Shoykhet B 1997. Micromechanics of diffusional creep. Int. J. Eng. Sci. 35, 1003–1032. (doi:10.1016/S0020-7225(97)00005-0)
8. Garikipati K, Bassman L, Deal M 2001. A lattice-based micromechanical continuum formulation for stress-driven mass transport in polycrystalline solids. J. Mech. Phys. Solids 49, 1209–1237. (doi:10.1016/S0022-5096(00)00081-8)
9. Villani A, Busso EP, Forest S 2015. Field theory and diffusion creep predictions in polycrystalline aggregates. Modelling Simul. Mater. Sci. Eng. 23, 055006 (doi:10.1088/0965-0393/23/5/055006)
10. Mishin Y, Waren JA, Sekerka RF, Boettinger WJ 2013. Irreversible thermodynamics of creep in crystalline solids. Phys. Rev. B 88, 184303 (doi:10.1103/PhysRevB.88.184303)
11. Gibbs JW. 1878. On the equilibrium of heterogenous substances. Trans. Conn. Acad. 3, 108–248, 343–524. Reprinted in The Scientific Papers of J. Willard Gibbs, 55–371. Oxbow Press.
12. Lowengrub J, Truskinovsky L 1998. Quasi-incompressible Cahn–Hilliard fluids and topological transition. Proc. R. Soc. Lond. A 454, 2617–2654. (doi:10.1098/rspa.1998.0273)
13. Lee EH. 1969. Elastic–plastic deformation at finite strains. J. Appl. Mech. 36, 1–6. (doi:10.1115/1.3564580)
14. Asaro RJ, Rice JR 1977. Strain localization in ductile single crystals. J. Mech. Phys. Solids 25, 309–338. (doi:10.1016/0022-5096(77)90001-1)
15. Mesarovic SDj, Padbidri J 2005. Minimal kinematic boundary conditions for simulations of disordered microstructures. Phil. Mag. 85, 65–78. (doi:10.1080/14786430412331313321)
16. Schneider T, Stoll E 1978. Molecular-dynamics study of a three-dimensional one-component model for distortive phase transitions. Phys. Rev. B 17, 1302–1322. (doi:10.1103/PhysRevB.17.1302)
17. Sutton AP, Balluffi RW 1995. Interfaces in crystalline materials. Oxford, UK: Clarendon Press.
18. Wolf D. 1989. Structure-energy correlation for grain boundaries in F.C.C. metals-I. Boundaries on the (111) and (100) planes. Acta Metall. 37, 1983–1993. (doi:10.1016/0001-6160(89)90082-5)
19. Haasen P. 1986. Physical metallurgy, 2nd edn Cambridge, UK: Cambridge University Press.
20. Souquet P. 1990. Une méthode simplifiée pour le calcul des propriétés élastiques de matériaux hétérogènes à structure périodique. C.R. Acad. Sci. 311, 769–774.
21. Mura T. 1991. Micromechanics of defects in solids, 2nd edn Dordrecht, the Netherlands: Kluwer Academic Publishers.
22. Moulinec H, Souquet P 1994. A fast numerical method for computing the linear and nonlinear properties of composites. C.R. Acad. Sci. 318, 1417–1423.
23. Carslaw HS, Jaeger JC 1959. Conduction of heat in solids, 2nd edn Oxford, UK: Oxford University Press.
24. Cammarata RC. 1994. Surface and interface stress effects in thin films. Prog. Surf. Sci. 46, 1–38. (doi:10.1016/0079-6816(94)90005-1)
25. Rice JR. 1981. Constraints on the diffusive cavitation of isolated grain boundary facets in creeping polycrystals. Acta Metall. 29, 675–681. (doi:10.1016/0001-6160(81)90150-4)
26. Anderson PM, Rice JR 1985. Constrained creep cavitation of grain boundary facets. Acta Metall. 33, 409–422. (doi:10.1016/0001-6160(85)90083-5)

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society