Lateral deformation of diffusion couples

Acta Materialia 53 (2005) 1995–2008 www.actamat-journals.com

Lateral deformation of diffusion couples W.J. Boettinger

a,*

, G.B. McFadden a, S.R. Coriell a, R.F. Sekerka

a,b

, J.A. Warren

a

a

b

National Institute of Standards and Technology, Gaithersburg, MD 20899, USA Departments of Physics and Mathematics, Carnegie Mellon University, Pittsburgh, PA 15213, USA Received 6 October 2004; received in revised form 7 January 2005; accepted 12 January 2005

Abstract A model is used to describe the shape change of a binary diffusion couple when the diffusivities of the two species differ. The classical uniaxial Kirkendall shift is obtained only if the displacement is constrained to be in the diffusion direction. For traction-free conditions at the external surfaces of a diffusion couple, a more general displacement field is obtained that accounts for the lateral shape change data of Voigt and Ruth [Journal of Physics-Condensed Matter 7 (1995) 2655–2666]. The model employs an isotropic stress-free strain rate and equal and constant partial molar volumes. In this case the displacement field is shown to be independent of the various elastic/plastic moduli. Depending on the lateral dimension of the diffusion couple, the displacement in the diffusion direction can be reduced by up to a factor of three compared to the case of a pure uniaxial displacement.
2005 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Bulk diffusion; Elastic behavior; Plastic deformation; Vacancies; Theory and modeling

1. Introduction The Kirkendall effect is a well-known consequence of a disparity in the intrinsic (lattice) diffusion coefficients of the components in alloy diffusion couples. The effect has three manifestations: a shift of inert markers parallel to the diffusion direction (deformation with pure axial displacement), lateral shape changes (displacements with a component orthogonal to the diffusion direction), and the formation of voids. These effects are summarized by Philibert [1] and may be understood for crystals in terms of the creation and elimination of vacancies in the diffusion zone. Current interest in stress generation in thin films due to diffusion processes has stimulated new interest in this classical problem. The usual 1-D analysis of interdiffusion occurring by a vacancy mechanism assumes that the material has a

*

Corresponding author. Tel.: +1 301 975 6160. E-mail address: [email protected] (W.J. Boettinger).

sufficient number of defects (dislocations and grain boundaries) to prevent a supersaturation of vacancies and to prevent void formation. This analysis leads to the prediction that lattice planes (with their associated inert markers) are displaced with respect to an observer embedded in the sample at a position outside of the diffusion zone. The velocity of the markers is parallel to the diffusion direction and is proportional to the vacancy flux with respect to the lattice. Philibert notes in his monograph that the lattice velocity in the diffusion (axial) direction may be a factor a (1/3 < a < 1) smaller due to the details of how the vacancies are annihilated and created. For example if vacancies are created/annihilated by the climbing of dislocations only in directions orthogonal to the diffusion direction (see Fig. 1), then, on average, the lattice velocity would be purely parallel to the diffusion direction and a would be equal to unity. On the other hand, if dislocations are present in a variety of orientations, then climb would occur in directions not necessarily orthogonal to the diffusion direction and a would not be unity. Clearly for an intrinsic dislocation

1359-6454/$30.00
2005 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2005.01.011

1996

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

r

z

Fig. 1. Schematic diagram of lattice planes and dislocations in a diffusion couple. Vacancies flow from right to left due to an assumed disparity of the atomic fluxes. The dislocations act as sinks (left) and sources (right) for the vacancies, as atoms are removed from the extra planes on the left and added to the planes on the right. The resultant climb of the dislocations causes, on average, a displacement of the lattice planes to the left indicated by the arrow marked with the letter V. This figure is adapted from Philibert [1].

density and distribution of Burger
s vectors typical of metals (and also grain orientations), a purely uniaxial displacement field is unlikely. In their early work on the Kirkendall effect, Correa da Silva and Mehl [2] were concerned with the possibility of lateral dimension changes. They concluded that the 1-D approximation appeared valid in their experiments. Huntington and Grones [3] studied electromigration in wires that were thin compared to the diffusion distance. They introduced the parameter a and assumed that dimensional change due to vacancy creation and annihilation would occur equally in all directions, making a = 1/3. Conversely, they argued that for samples with thicknesses that are large compared to the diffusion distance, like those of Correa da Silva and Mehl, the shape change would occur mostly in the diffusion direction, making a nearly unity. The present paper recovers these limiting cases. Fig. 2 shows a drawing of a classical unidirectional Kirkendall displacement as well as a hypothetical displacement field with lateral shape change. The outward bulge would occur on the side of the diffusion couple rich in the slow diffusing species. Ruth and coworkers [4–7] have performed detailed measurements of the lateral shape of a series of Ag–Au diffusion couples. They [6,7] have modeled the distortion by assuming that a is unity from the centerline of the sample out to a critical distance from the lateral surface where it changes to one half (in a slab geometry). Estimates for this critical distance are based on the shear strength of the material. In this paper, Section 2 presents the governing equations for a 3-D model that recovers the 1-D Kirkendall

Fig. 2. Drawing depicting two types of deformations that might occur in a diffusion couple depending on how vacancies are annihilated or generated. The diffusion direction is along the z axis and vacancies are generated on the left and annihilated on the right. The top figure shows a system of markers prior to diffusion. The middle figure shows the classical unidirectional Kirkendall displacement field. The bottom figure shows a hypothetical deformation that approaches the unidirectional Kirkendall behavior in the interior but differs toward the free surfaces. This figure is adapted from one shown by Philibert [1].

displacement field as a special case, develops a solution method using normal modes for the deformation of cylindrical and slab diffusion couples with zero-traction boundary conditions on the lateral surfaces, and gives limiting behavior for thin and thick cylinders and slabs. Section 3 gives numerical results of the model for the lateral shape change which compare favorably with the experimental results of Voigt and Ruth. Lateral bulging is obtained naturally from the model with no assumptions about the value of a. Numerical results are also presented that describe the displacement field throughout the diffusion couple.

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

1997

2. Model

neglected, we recover the equation of mechanical equilibrium,

2.1. General 3-D model

orjk ¼ 0; oxk

We seek a model that does not treat the discreteness of the lattice or a detailed spatial description of events around defects; e.g., a climbing dislocation. Thus we consider a model that employs a continuum, with field variables that average over volumes large enough to contain many randomly oriented defects, but small compared to diffusion distances. As such, we employ a continuum description similar to that used by Stephenson [8]. Positions in the initially-undeformed reference state are denoted by x, and a displacement field u(x,t) is defined so that a coarse-grained lattice point x at time t = 0 is at position x + u(x,t) at time t. In the small strain approximation, the total strain tensor eTij , given by
 1 oui ouj eTij ¼ þ ð1Þ 2 oxj oxi is separated into three parts: stress-free, viscous (plastic), and elastic, with eTij ¼ e0ij þ ePij þ eEij . The time derivative of the displacement field gives the Lagrangian flow velocity v(x,t) = ou/ot, and the rate of strain tensor is then
 1 ovi ovj þ ; ð2Þ e_ Tij ¼ 2 oxj oxi where e_ Tij ¼ oeTij =ot denotes the time derivative of the total strain. The stress-free strain rate e_ 0ij is related to the vacancy creation rate rv and the partial molar volume of vacancies, V v . Following the coarse graining assumption, the stress-free rate of strain tensor, i.e., the rate of strain expected with no stress and/or geometrical constraint of the sample is assumed to be isotropic and is given by [8] e_ 0ij ¼ 13dij V v rv :

ð3Þ

As will be seen, the total strain will not be isotropic. Assuming a Newtonian viscous material, the plastic strain rate is a function of the stress by the relationship   1 1 rij
dij rkk ; e_ Pij ¼ ð4Þ 2g 3 where g is the viscosity. Here, we employ a summation convention on repeated indices. The plastic flow is thus assumed incompressible, as e_ Pkk ¼ 0. The isotropic elastic strain is given by eEij ¼

1þm m rij
dij rkk ; E E

ð5Þ

where m is Poisson
s ratio and E is Young
s modulus. To simplify the computation, we will assume that m = 1/2, so that the elastic strain is also incompressible, eEkk ¼ 0. If we assume inertial effects and body forces can be

j ¼ 1; 2; 3:

ð6Þ

For closure we employ standard equations for the vacancy source (see, e.g. [1, p. 307]), V v rv ¼
r  ½ðDA
DB ÞrX B 

ð7Þ

and the diffusion field (see, e.g. [1, p. 209]), oX B ¼ r  ½ðDA X B þ DB ð1
X B ÞÞrX B ; ot

ð8Þ

where DA and DB are the intrinsic diffusion coefficients of components A and B, and XB is the mole fraction of component B. A crucial simplification in this work is that we have assumed that the partial molar volumes of components A, B, and vacancies are identical and constant; specifically, the partial molar volumes are independent of stress, as reflected in our choice of m = 1/2. With the assumption of equal partial molar volumes, the atomic fluxes are independent of stress. Moreover, in writing Eq. (8) we have omitted nonlinear terms that would couple the concentration gradient to the deformation field. These assumptions give rise to a Kirkendall effect that results only from the difference of intrinsic diffusivities, and lead to a more tractable model analytically. 2.2. Solution procedure for traction-free problems We next outline a solution procedure for calculating the displacement fields for cylindrical and slab geometries with zero traction boundary conditions for a given solution to the diffusion problem. A combination of definitions and constitutive relations can be written as

 1 ovj ovk 1 o oX B þ D0 þ djk 2 oxk oxj 3 oxp oxp

 1 1 3 1 rjk
rpp djk þ ¼ r_ jk
r_ pp djk ð9Þ 2g 3 2E 3 for j, k = 1, 2, 3, where D0 = DA
DB. We also need to specify appropriate initial and boundary conditions on the system. The trace of the right hand side of Eq. (9) vanishes, so that we have r  v þ r  ½D0 rX B  ¼ 0;

ð10Þ

where $ Æ v is the divergence of the velocity field. If we let

 1 3g 2g o oX B rpp þ r_ pp
P¼
D0 ð11Þ 3 E 3 oxp oxp denote an effective pressure field, we may write Eq. (9) in the form

1998

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008


 3g ovj ovk rjk þ r_ jk ¼
P djk þ g þ ; E oxk oxj

j; k ¼ 1; 2; 3: ð12Þ

By differentiating Eq. (6) with respect to time we also have or_ jk ¼ 0; oxk

j ¼ 1; 2; 3

ð13Þ

and so orjk 3g or_ jk þ ¼ 0; E oxk oxk

j ¼ 1; 2; 3:

ð14Þ

Using Eqs. (14) and (12), we therefore obtain the equation rP ¼ gr2 v þ grðr  vÞ;

ð15Þ

which is a (compressible) version of the Stokes equation from hydrodynamics. This suggests the use of solution techniques that have been well-developed in the hydrodynamics literature [9]. By using the identity $($ Æ v) = $2v + $ · ($ · v), Eq. (15) can be written as rP ¼
gr  ðr  vÞ þ 2grðr  vÞ:

ð16Þ

The velocity field satisfying Eq. (10) can be expressed in the form v ¼
D0 rX B
r  a;

ð17Þ

where a is a vector potential. Taking the divergence of the velocity recovers Eq. (10), and taking the curl gives r  v ¼
r  ðr  aÞ:

ð18Þ

Here, we have used the fact that $D0 · $XB vanishes since D0 is only a function of XB. Taking the divergence and curl of Eq. (16) then yields r2 P ¼ 2gr2 ðr  vÞ ¼
2gr2 ðr  ½D0 rX B Þ

ð19Þ

and r  ½r  fr  ðr  aÞg ¼ 0:

ð20Þ

Given a diffusion solution XB, Eqs. (19) and (20) determine P and a that are then used to calculate the velocity and displacement fields. Boundary conditions for these equations are obtained from the condition of vanishing surface traction, rjknk = 0. Consequently we also have r_ jk nk ¼ 0 at the surface, and so   3g rjk þ r_ jk nk ¼ 0; ð21Þ E where we have employed n_ k ¼ 0, since the surface in the reference frame is stationary. By using Eq. (12), the traction boundary condition can be expressed in terms of P and $v. It is notable that the solution for the velocity and displacement fields does not depend on the values of E and g; however, the associated stress field does depend on these parameters.

2.3. One-dimensional diffusion model In order to obtain a tractable problem, we assume that the diffusion involves sufficiently small concentration dife ¼ DA X B þ DB ð1
X B Þ ferences that the parameters D and D0 = DA
DB that appear in Eqs. (7) and (8) are constants. A concentration profile XB(z,t) representing a one-dimensional diffusion couple for
1 < z < 1 then has the form ! L L ðX R ðX R z B þ X BÞ B
X BÞ þ erf pffiffiffiffiffiffiffiffi ; X B ðz; tÞ ¼ ð22Þ 2 2 e 4 Dt where X LB and X R B are the far field atomic fractions, and we have written z = x3 to simplify the notation. Here Z z 2 2 p ffiffiffi erfðzÞ ¼ e
s ds ð23Þ p 0 is the error function. The stress-free strain rate tensor is then given by D 0 o2 X B djk 3 oz2
 L 2 D0 ðX R z2 B
X BÞ ¼ pffiffiffi e 3=2 z exp
e djk : 3 4 Dt pð4 DtÞ

e_ 0jk ¼


ð24Þ

This stress-free strain rate for the above one-dimensional diffusion field will be used for all subsequent calculations. 2.4. Displacement fields We consider analytical solutions for the velocity and displacement fields, given the time-dependent concentration profile XB(z,t) and the associated stress-free strain rate. We start with a solution with a purely axial displacement field (not traction-free) that recovers the classical Kirkendall effect. We then treat higherdimensional solutions in traction-free geometries that allow displacements normal to the diffusion gradient. We consider a cylindrical geometry, corresponding to a reference state given by an infinite rod of radius d, with
1 < z < 1 and 0 < r < d. For comparison purposes we then consider a slab of width 2d, with the reference state
d < y < d and
1 < z < 1; the displacement is assumed to vanish in the third dimension (plane strain). 2.4.1. Pure axial displacement For a single non-zero displacement uz(z,t) driven by the diffusion profile XB(z,t) given in Eq. (22), the velocity vz(z,t) = ouz(z,t)/ot must satisfy Eq. (10), ovz o2 X B ¼
D0 ; oz oz2

ð25Þ

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

which gives e

ð26Þ

:

The corresponding axial displacement is given by pffiffiffiffiffiffiffiffiD  ðX R
X L Þ 0 B B e uz ðz; tÞ ¼ 4 Dt e 2 D ! " # jzj jzj 1
z2 =4e Dt  pffiffiffiffiffiffiffiffi erfc pffiffiffiffiffiffiffiffi
pffiffiffi e ; ð27Þ p e e 4 Dt 4 Dt where erfc(z) = 1
erf(z). This means that markers initially at position z have, after a time t, moved to the position z + uz(z,t). The Kirkendall plane is therefore located at the point uz(0,t), which is given by pffiffiffiffiffiffiffiffiD  ðX R
X L Þ 0 B e pffiffiffi B ; uz ð0; tÞ ¼
4 Dt ð28Þ e 2 p D pffi so that the displacement at z = 0 goes as t. At finite times, the displacement for z > 0 is monotonic in p z ffiffiffiffiffiffiffiffi and e = 4 Dt vanishes at infinity. The scaled displacement u pffiffiffiffiffiffiffiffi z e . These reis a function of the similarity variable z= 4 Dt sults are all well-known. For future reference, we note that the concentration field in Eq. (22) and the displacement field in Eq. (27) can be written in Fourier space as [10] L ðX R ðX R
X LB Þ B þ X BÞ þ B 2 p

Z 0

1

sin kz e expð
k 2 DtÞdk; k ð29Þ

  R D0 ðX B
X LB Þ uz ðz; tÞ ¼
e p D Z 1 cos kz e  f1
expð
k 2 DtÞgdk: 2 k 0

6

5

4

3

2

1 0

0.2

0.4

0.6

0.8

1

XAg

Fig. 3. The solid curve shows the interdiffusion coefficient, e ¼ DAu X Ag þ DAg ð1
X Ag Þ, and the dashed curve shows D
D0 = DAg
DAu, as functions of the mole fraction of silver, XAg. DAg and DAu are determined using [11–13].

In Fig. 4, we show plots of XB(z,t) and V v rv , together with the one-dimensional velocity and axial displacement fields using these diffusion coefficients. The plots are shown for the maximum experimental time used by Ruth et al. of 8100 min (t = 4.86 · 105 s).

0.3

XB

X B ðz; tÞ ¼

7

0.2 0.1 0 –1

ð30Þ

As we are ultimately interested in computing lateral deformation corresponding to the experiments of Ruth et al., we require values for D0 and the interdiffusion e for Au–Ag alloys at 1081 K. Fig. 3 shows coefficient D the values obtained from assessed thermodynamic [11] and diffusion mobility [12] databases via the DICTRA code [13]. With XB = XAg, we evaluate the diffusion coefficients for a concentration at the midpoint (XAg = 1/8) of one of the diffusion couples considered by Ruth et al. with pure Au against Au–25% Ag. From Fig. 3, we find that D0 ¼ DAu
DAg ¼
3:3  10
14 m2 =s;

ð31Þ

e ¼ DAu X Ag þ DAg ð1
X Ag Þ ¼ 7:2  10
14 m2 =s: D

ð32Þ

We note that our assumption that the partial molar volumes of the two species are equal is true to a good approximation for the Ag–Au system, for which the molar volumes of pure Ag and pure Au differ by about one percent.

108 Vv σv (s–1 )


pffiffiffi p

105 vz (µm/s)

X LB Þ
z2 =4e Dt

1014 Diffusion Coefficients (m/s2)

8

ðX R B

–0.5

0

0.5

1

–0.5

0

0.5

1

–0.5

0

0.5

1

0

0.5

1

4 2 0 –2 –4 –1 1.5 1 0.5 0 –1 16

uz (µm)


D0 vz ðz; tÞ ¼ pffiffiffiffiffiffiffiffi e 4 Dt

1999

12 8 4 0 –1

–0.5

z (mm) Fig. 4. Profiles of the error function concentration profile XB for e ¼ 7:2  10
14 m2/s, the vacancy creation/annihilation rate V v rv for D D0 =
3.3 · 10
14 m2/s, the one-dimensional velocity vz and displacement uz for t = 4.86 · 105 s.

2000

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

We note that the discontinuity in slope of uz at z = 0 in Fig. 4 is an artifact due to the discontinuity of XB at z = 0 and t = 0 of the diffusion couple as given by Eq. (22). This causes vz in Eq. (26) to be proportional to a delta function of z at t = 0. One can eliminate this discontinuity by replacing vz(z,t) in Eq. (26) by vz(z,t + t0), which is equivalent to starting the calculation after a small amount of diffusion has taken place, which experimentally will always occur during heating. Eq. (27) is then replaced by the function uz(z,t + t0)
uz(z,t0), which has no discontinuity of slope at any time.

The normal to the surface r = d is ^r, and the boundary conditions (21) with Eq. (12) take the form

2.4.2. Displacement field in a cylindrical geometry In cylindrical coordinates (r, h, z), with z = x3, we look for an axisymmetric solution by writing the vector potential in terms of a Stokes stream function w(r, z), with

ð42Þ

1 h: a ¼ wðr; zÞ^ r

ð33Þ

Eq. (17) then takes the form v¼

1 ow 1 ow oX B

D0 : r oz r or oz

ð34Þ

An equation for the stream function follows from Eq. (20) in the form 

2 o2 1 o o2 þ
w ¼ 0: or2 r or oz2

½rbv rz þ ½rbv zr ¼ 0;

ð41Þ

where $v is the tensor given by the gradient of the deformation rate in cylindrical coordinates. Applying the boundary conditions, we obtain A¼

B¼


kdD0 X 0 ½I 0 ðkdÞ þ I 2 ðkdÞ ½I 0 ðkdÞ þ I 2 ðkdÞ½I 1 ðkdÞ þ kdI 0 ðkdÞ
2kd½I 1 ðkdÞ

1

X 0 ðkÞe


k 2 e Dt

sin kz dk;

: 2 ½I 0 ðkdÞ þ I 2 ðkdÞ½I 1 ðkdÞ þ kdI 0 ðkdÞ
2kd½I 1 ðkdÞ ð43Þ

L ðX R B
X BÞ : ð37Þ pk We may then construct solutions by a normal mode analysis; the stream function has modes proportional to cos kz, and the effective pressure has modes that are proportional to sin kz. The full solution in physical space is then given by integration of the modes over the wavenumber k for 0 < k < 1. The normal modes are denoted by circumflexes. The normal mode solution to Eq. (35) for the stream function is given by Dt b zÞ ¼ e
k2 e wðr; r½AI ðkrÞ þ BkrI ðkrÞ cos kz; ð38Þ 2

where the constants A and B are determined by imposing the zero traction boundary conditions at r = d, and In(kr) denotes the modified Bessel function of the first kind of order n [14]. From Eq. (16) the corresponding effective pressure mode is found to have the form Dt b zÞ ¼ 2gk 2 e
k2 e Pðr; ½D0 X 0 ðkÞ
BI 0 ðkrÞ sin kz:

ð39Þ

ð44Þ

2 Dt ^vz ¼
ke
k e ½D0 X 0 ðkÞ þ AI 0 ðkrÞ þ krBI 1 ðkrÞ cos kz:

ð45Þ Performing a time integration, we obtain the normal modes for the displacement field

^ uz ¼

X 0 ðkÞ ¼

;

2D0 X 0 I 1 ðkdÞ

2 Dt ^vr ¼
ke
k e ½AI 1 ðkrÞ þ krBI 2 ðkrÞ sin kz;

ð36Þ

0

2

In these expressions the denominator is positive and monotonically increasing in k. The resulting normal modes for the velocity field have the form

^ur ¼

Z

with

1

ð40Þ

ð35Þ

Our approach is to use a Fourier representation in z, in which case Eq. (29) leads to ðX R þ X LB Þ þ X B ðz; tÞ ¼ B 2

b þ 2g½rbv  ¼ 0;
P rr


1 2 Dt g½AI 1 ðkrÞ þ krBI 2 ðkrÞ sin kz; f1
e
k e e kD

ð46Þ


1 2 Dt g½D0 X 0 ðkÞ þ AI 0 ðkrÞ þ krBI 1 ðkrÞ cos kz; f1
e
k e e kD ð47Þ

where we have chosen the integration constant to be consistent with zero initial displacement. Evaluating along the axis r = 0 gives ^ur ¼ 0, and ^uz ¼


L ðX R B
X B Þ D0 cos kz e f1
expð
k 2 DtÞgh z ðkdÞ; e k2 p D ð48Þ

where hz ðkdÞ " ¼ 1


#

kd½I 0 ðkdÞ þ I 2 ðkdÞ ½I 0 ðkdÞ þ I 2 ðkdÞ½I 1 ðkdÞ þ kdI 0 ðkdÞ
2kd½I 1 ðkdÞ2

:

ð49Þ

Evaluating at the surface r = d gives  L  ðX R 2 Dt B
X B Þ D0 sin kz ^ur ¼ f1
e
k e ggr ðkdÞ; 2 e p k D

ð50Þ

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

where gr ðkdÞ ¼

2½I 1 ðkdÞ

½I 0 ðkdÞ þ I 2 ðkdÞ½I 1 ðkdÞ þ kdI 0 ðkdÞ
2kd½I 1 ðkdÞ2 ð51Þ

and ^ uz ¼


2

ðX R B


p

X LB Þ

  D0 cos kz 2 Dt f1
e
k e ggz ðkdÞ; 2 e k D

ð52Þ

where gz ðkdÞ ¼

I 1 ðkdÞ½I 0 ðkdÞ þ I 2 ðkdÞ

: ½I 0 ðkdÞ þ I 2 ðkdÞ½I 1 ðkdÞ þ kdI 0 ðkdÞ
2kd½I 1 ðkdÞ2 ð53Þ

To obtain the components of the displacement field in physical space we then integrate the modes ^ ur and ^uz with respect to the wavenumber k. 2.4.3. Limiting cases of thick or thin samples (cylindrical case) pffiffiffiffiffiffi e and the dimensionless The diffusion length scale Dt e 2 play fundamental time Dt=d roles in this process. A pffiffiffiffiffiffi e Þ then correspondspto thick sample ðd  Dt ffiffiffiffiffiffi short e Þ cordimensionless times, and a thin sample ðd  Dt responds to long dimensionless times. For discussion purposes p weffiffiffiffiffirefer to thick ffi pffiffiffiffiffiffior thin samples as meaning e or d  Dt e , respectively. that d  Dt For thick samples the dominant contribution to the Fourier integrals comes from the region kd  1. In Eqs. (48)–(52) this limiting case corresponds to setting hz(kd), gr(kd), and gz(kd) to unity. The resulting forms for the axial normal mode are particularly simple for either r = 0 or r = d, and are both given by the normal mode in the onedimensional Fourier transform result (30). Consequently, the corresponding axial component is the same as the one-dimensional case given by Eq. (27) for both the displacement at r = 0 and r = d. The limiting form of the radial normal mode at r = d does not have a simple closed-form solution in physical space, and we have resorted to numerical integration as discussed in Section 3. The displacements for thin samples correspond to the limit kd  1. By expanding the Bessel functions in Eqs. (48)–(53) for small kd, the thin limit of the normal modes in this case gives, to leading order, ^ uz ¼


L ðX R 2 Dt B
X B Þ D0 cos kz f1
e
k e g e k2 3p D

ð54Þ

for both r = 0 and r = d; as shown below, this equation holds for all r. This displacement is 1/3 of the displacement for the thick limit and the one-dimensional case, a result consistent with the argument of Huntington and Grones [3]. The limiting form of the radial normal mode at r = d is, to leading order, ^ ur ¼ d

L ðX R 2 Dt B
X B Þ D0 sin kz f1
e
k e g: e k 3p D

ð55Þ

2001

In this case the integration can be performed in closed form and gives the radial displacement field of the outer surface !  L  ðX R jzj B
X B Þ D0 erfc pffiffiffiffiffiffiffiffi sgnðzÞ; ur ðz; d; tÞ ¼ d e 6 e D 4 Dt

ð56Þ

where sgn(z) = z/jzj is the sign of z for z 6¼ 0. This limiting solution is discontinuous at z = 0, and will be investigated in more detail below. Alternatively, the displacements in the thin limit can be obtained directly from the governing equations using scaling arguments similar to those in lubrication theory or shell theory (see, e.g., [15]). To simplify the discussion of the boundary conditions relevant to the thin limit, we consider the case of elasticity without plastic flow. We emphasize that the same results hold for the more general case if the stress tensor rjk is replaced by the tensor rjk þ ð3g=EÞr_ jk as described in Section 2.2. For a thin cylindrical sample, the traction-free boundary conditions rrr = 0 and rrz = 0 at r = d are effectively impressed throughout the entire sample. The remaining equilibrium condition, orzz/oz = 0, can be integrated directly and, since the sample is stress-free in the far-field where z ! ±1, the entire sample is stress-free in the thin limit. Inffi the thin limit the appropriate scale for vz and z is pffiffiffiffiffi e , while both vr and r scale with d. Thus, we have Dt jovr/ozj  jovz/orj, and the equation rrz = 0 implies that to leading order vz = vz(z,t) is independent of r. From rzz = 0 it follows that the effective pressure is also independent of r and is given by Pðz; tÞ ¼ 2govz =oz. The equation rrr = 0 then gives vr ¼ rPðz; tÞ=2g. From Eq. (10), we then find that 2 o2 X B ; Pðz; tÞ ¼
gD0 3 oz2

ð57Þ

and so we have vz ¼


D0 oX B ; 3 oz

vr ¼


rD0 o2 X B : 3 oz2

ð58Þ

This is a general result for any profile XB(z,t). Inserting the expression for XB from Eq. (22) and performing a time integration to obtain the displacement fields from the velocity fields gives pffiffiffiffiffiffiffiffi   e D0 ðX R
X L Þ 4 Dt B B uz ðz; tÞ ¼ e 2 3 D ! " # jzj jzj 1
z2 =4e Dt  pffiffiffiffiffiffiffiffi erfc pffiffiffiffiffiffiffiffi
pffiffiffi e ; p e e 4 Dt 4 Dt

ð59Þ

2002

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

!  L  ðX R jzj B
X B Þ D0 erfc ur ðz; r; tÞ ¼ r sgnðzÞ: e e 1=2 6 D ½4 Dt ð60Þ The latter expression reduces to Eq. (56) for r = d, and the former expression is again 1/3 of the one-dimensional result. We note that these limiting expressions hold throughout the thin sample. 2.4.4. Displacement field in a slab geometry We consider a two-dimensional solution for a slab geometry given by
1 < z < 1 and
d < y < d. We assume there is no displacement in the x direction, and that there are traction-free boundary conditions at y = ±d. We again represent the flow field in the form of Eq. (17), where we now write a ¼ wðy; zÞ^x:

ð61Þ

This gives vy ¼


ow ; oz

vz ¼

ow oX B
D0 : oy oz

b z; tÞ ¼ cos kz½E sinh ky þ Fky cosh ky expð
k 2 DtÞ; e wðy; ð63Þ where E and F are coefficients that are obtained by applying the traction-free boundary conditions at y = d. The corresponding effective pressure mode has the form 2

b z; tÞ ¼ 2gFk 2 sin kz cosh ky expð
k 2 DtÞ e
2gD0 o X B : Pðy; oz2 ð64Þ Applying the no-stress boundary conditions at y = d gives D0 X 0 ½sinh kd þ kd cosh kd ; ½cosh kd sinh kd þ kd


D0 X 0 sinh kd F ¼ : ½cosh kd sinh kd þ kd

^ uz ¼

D0 cos kz e f1
expð
k 2 DtÞg e pk 2 D sinh kd cosh kd ;  ½cosh kd sinh kd þ kd

L ^uz ¼
ðX R B
X BÞ

D0 sin kz e f1
expð
k 2 DtÞg e pk 2 D sinh2 kd :  ½cosh kd sinh kd þ kd

ð70Þ

L ^uy ¼ ðX R B
X BÞ

ð71Þ

2.4.5. Limiting cases of thick or thin samples (slab geometry) As in the cylindrical geometry, the dependence of the solutions on the sample thickness can be eliminated by choosing dimensionless variables based on the length e We may again consider scale d and the time scale d2 = D. the cases of p thick and thin samples relative to the diffuffiffiffiffiffiffi e . The thick limit kd  1 gives the same sion length Dt limiting forms for the normal modes as in the cylindrical case, where the mode ^uy now plays the role of the radial mode ^ur . The displacement fields for thick samples at the outer surface and at the mid-section are therefore independent of the specific geometry. In contrast, the thin limit of the normal modes is obtained by expanding the hyperbolic functions in Eqs. (69)–(71) for small kd, giving, to leading order,

ð65Þ

ð66Þ

1 e sinh ky þ Fky cosh ky; sin kzf1
expð
k 2 DtÞg½E e kD ð67Þ 1 e cos kzf1
expð
k 2 DtÞg e kD  ½ðE þ F Þ cosh ky þ Fky sinh ky
D0 X 0 ;

ð69Þ

At y = d we have

^uz ¼


The resulting normal modes for the displacement are ^ uy ¼

D0 cos kz e f1
expð
k 2 DtÞg e pk 2 D ½sinh kd cosh kd þ kd
kd cosh kd :  ½cosh kd sinh kd þ kd

L ^uz ¼
ðX R B
X BÞ

ð62Þ

The appropriate periodic normal mode, which is odd about the plane y = 0, has the form

E¼

where we have chosen the integration constant to be consistent with zero initial displacement. These expressions simplify somewhat at the centerline, y = 0, and at the boundary, y = d. At y = 0 we have ^uy ¼ 0 and

ð68Þ

L ðX R 2 Dt B
X B Þ D0 cos kz f1
e
k e g e k2 2p D

ð72Þ

for both y = 0 and y = d, which is 1/2 of the result for the thick limit. The limiting form of the normal mode uy at y = d is, to leading order, ^uy ¼ d

L ðX R 2 Dt B
X B Þ D0 sin kz f1
e
k e g: e 2p k D

ð73Þ

This displacement is of order d, and vanishes in the thin sample limit. Unlike the case of thick samples, for thin samples the sample geometry affects the limiting forms of the displacement at the outer surface and at the mid-section. The axial displacement in the cylindrical geometry is a factor of 1/3 of the axial displacement for a thick sam-

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

2003

ple, whereas for the slab geometry it is a factor of 1/2 of the thick sample
s displacement. As in the cylindrical case, it is also possible to obtain the displacements in the thin limit directly from the governing equations using thin plate theory. The results show the same geometrical factor of 1/2, rather than 1/3, appearing in the resulting displacement fields; the effective pressure field has the same form in either geometry.

3. Numerical results for the cylindrical case The displacements were obtained in the quadrant z > 0 and r > 0 by numerical calculation of the Fourier transforms using the QUADPACK software QAWF by Piessens et al. [16]. The z-displacements are symmetric about r = 0 and about z = 0. The radial displacements are symmetric about r = 0 and antisymmetric about z = 0. 3.1. Dimensional results and comparison to experiment Fig. 5 shows plots of the radial and axial displacement fields for all four quadrants for a 0.5 cm radius cylindrical diffusion couple for the same parameters as were used in Fig. 4. Note that the axial displacement as a function of z has the same peaked shape as the 1-D solution except near the cylinder surface. The radial displacement is small except near the cylinder surface where it is positive and negative on opposite sides of z = 0. This radial displacement produces the change in shape of the exterior surface. Thus along the center section of the rod the displacement field is similar to the 1-D result. Ruth et al. have measured the external surface shape for a diffusion couple with X LAg ¼ 0 and X R Ag ¼ 0:25 annealed at 1081 K for 4.86 · 105 s. The specimens were half cylinders with semi-circular cross-sections that are 0.5 cm in radius. As shown in Fig. 6, the measured profile is slightly asymmetric, with a bulge of about 5.0 lm, and a depression of about 6.3 lm, with a peak-to-peak separation of about 265 lm. We overlay the surface shape obtained from the calculations shown in Fig. 5. The displacement field gives a deflection of ur = 5.3 lm for both the bulge and the depression, with a peak-topeak separation of 336 lm; the maximum value of ouz/or is approximately 0.1. The calculated profile is the curve given by the Cartesian points [d + ur(d,z,t),z + uz(d,z,t)] as z ranges over the reference configuration. Because Ruth et al. did not precisely identify the location of the initial diffusion interface, we have shifted the computed curves in z so that the curves coincide at the points where the surface defore 2 ¼ 1:4 mation The dimensionless time Dt=d pvanishes. ffiffiffiffiffiffi
3
2 e =d ¼ 3:7  10  for the experiment is small 10 ½ Dt

Fig. 5. Computed axial, uz, and radial, ur, displacements versus position (r,z) in a 1 cm diameter cylindrical diffusion couple at e ¼ 7:2  10
14 m2/s. t = 4.86 · 105 s for D0 =
3.3 · 10
14 m2/s and D

Fig. 6. A comparison of the measured bulge and the calculated bulge, ur(d,z) [dashed], and ur(d,z + uz) (solid).

enough that the computed curve for the cylinder with d = 0.5 cm is graphically indistinguishable from the limiting curve for cylindrical samples of infinite radius, and is also indistinguishable from the curve for a slab

2004

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

geometry with d = 0.5 cm. Thus the semi-circular geometry of the experiment is not relevant. The agreement between the measured and calculated lateral shape change of the diffusion couple is a major accomplishment of this paper. To our knowledge other aspects of the deformation away from the centerline of diffusion couples have not been measured. In Fig. 7, the computed axial displacement uz versus r at z = 0 and at z = 168 lm is compared to the 1-D result. The computed displacements are slightly smaller than the 1-D result at both values of z. These parameter values correspond to the thick limit, and, as discussed in Section 2.4.1, the axial displacements are approximately the same at r = 0 and r = d. For intermediate values of r the axial displacement reaches a minimum value which is localized near the surface r = d. In the thick limit, the pffiffiffiffiffi ffi only remaining length e . For t = 4.86 · 105 s, scale Dt p ffiffiffiffiffiffi in the problem is e  0:02 cm is roughly the distance from the surface Dt to the minimum value. Thus, the simple nature of ur versus r shown in the lower part of the schematic in Fig. 2 is not correct in detail. The shape of the ur versus r curve for fixed z can be even more complex than appears in Fig. 7. The condition (41) of vanishing tangential stress at the surface implies that our/oz + ouz/or = 0 at r = d. Thus at the extremum in ur along the surface the quantity ouz/or undergoes a sign change. This can be seen from careful examination of the plot of uz in Fig. 5 near r = d. This inflection is difficult to observe in Fig. 7 due to the scale, but has been verified by more detailed calculations near r = d. Whether these particulars indicate an observable behavior or are merely a result of the model assumptions is not known. The corresponding radial displacement field ur for 0 < r < d near the bulge minimum at z = 168 lm is much

simpler, as shown Fig. 8. The magnitude of the radial displacement increases monotonically from zero at r = 0 to 5.3 lm at r = d. The radial displacement field at z = 0 is zero for all r. 3.2. Comparison of numerical results in dimensionless form to thick and thin limiting cases If the coordinates and displacement fields are scaled e the resulting dimensionless with d, and time with d2 = D, solution is independent of the sample thickness; viz., we can define dimensionless displacement fields Ur and Uz in the cylindrical case for z > 0 for various dimensionless e 2 via times t D=d  L
ur ðr; z; tÞ ðX R B
X B Þ jD0 j e 2 Þ; ð74Þ ¼ U r ðr=d; z=d; t D=d e d p D  L
uz ðr; z; tÞ ðX R B
X B Þ jD0 j e 2 Þ: ð75Þ ¼ U z ðr=d; z=d; t D=d e d p D We then look at the behavior for small and large vale 2. ues of t D=d 3.2.1. Thin limit As shown in Fig. 9 the axial displacement at r = d (solid curves) increases steadily with time. The displacement is largest at z = 0, and decays for large z with a length scale that increases in time. We also show the axial displacement at r = 0 (dotted curves). The displacement differs slightly from that for r = d. For the range of times considered in this figure, the 2-D solution is described by the thin limit approximation, and the axial displacements are about one-third the size of the onedimensional axial displacements (dashed curves). The approximation is worst near z = 0, where there is a sys-

0.5

0.5 0.4

0.3

r (cm)

r (cm)

0.4

0.2

0.2

0.1 0 2

0.3

0.1

4

6

8 uz (µm)

10

12

14

Fig. 7. The axial displacement uz for 0 < r < d at z = 0 (solid curve) and near the bulge minimum at z = 168 lm (dashed curve) for d = 0.5 cm at t = 4.86 · 105 s. The corresponding one-dimensional displacements are shown as dotted lines. The independent variable is plotted on the vertical axis for comparison to the schematic in Fig. 2.

0 –6

–5

–4

–3 ur (µm)

–2

–1

0

Fig. 8. The radial displacement ur for 0 < r < d near the bulge minimum at z = 168 lm (dashed curve) for d = 0.5 cm at t = 4.86 · 105 s; the radial displacement vanishes for all radial positions r at z = 0.

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008 7

0.1

6

0

5

2005

–0.1

4

Uz

–0.2 3

–0.3 2

–0.4

1

–0.5

0 –1 0

0.5

1

1.5

2

2.5

3

0

0.5

1

1.5

2

2.5

3

z/δ Fig. 9. The dimensionless axial displacement Uz for 0 < z/d < 3 for e 2 ¼ 0:1, 1.0, 10.0, and dimensionless times (from bottom to top) Dt=d 100.0. The solid curves show Uz at the outer surface r = d, the dotted curves show Uz at the centerline r = 0, and the dashed curves show the e 2  1, which is one-third corresponding thin limit approximation Dt=d of the one-dimensional axial displacement.

Fig. 10. The dimensionless radial displacement Ur at the outer surface for 0 < z/d < 3 for dimensionless times (from top to bottom) e 2 ¼ 0:1, 1.0, 10.0, 100.0, and 10,000. The solid curves indicate Dt=d the numerical calculation of the general solution, and the dashed curves, which issue from the value
p/6 at z = 0, show the corree 2  1. sponding approximate solutions in the thin limit, Dt=d

tematic deviation for both the surface displacement at r = d and the displacement along the centerline at e 2 has the limr = 0. The solution for large values of Dt=d iting behavior " # ! Uz p jzj jzj 1
z2 =4e Dt pffiffiffiffiffiffi ¼
pffiffiffiffiffiffiffiffi erfc pffiffiffiffiffiffiffiffi
pffiffiffi e ; 3 p e e e =d 4 Dt 4 Dt Dt ð76Þ pffiffiffiffiffiffi e where the right hand side depends only on z= Dt. The radial displacement of the surface shown in Fig. 10 exhibits a depression which deepens and moves outward from z = 0 as the diffusion takes place, and, at any fixed position, approaches a limiting finite displacement for long times. The approximate solutions for Ur at r = d defined by Eq. (56) in the thin limit are shown by the dashed curves in Fig. 10. The approximate solutions all share the value
p/6 at z p = ffiffiffiffiffi 0, ffiand as seen in Eq. (56) depend only on e . The limiting solution is discontinuthe variable z= Dt ous at z = 0; this behavior can be seen clearly in Fig. 11, where pffiffiffiffiffiffi the solutions are replotted versus the variable e . The solution curves for all finite d are, however, z= Dt continuous across z = 0. Note Ur has a limiting pffiffiffiffiffithat ffi e =d, whereas the ampliamplitude for large values of p Dt ffiffiffiffiffiffi e =d. tude of Uz grows linearly with Dt 3.2.2. Thick limit A close comparison of the computed axial and radial displacement fields with the small time (large d) approximate solution at t = 4.86 · 105 s shows that the approximation is less accurate for the axial solution than for the radial solution at this time. In Fig. 12, we show the dimensionless axial displacement Uz at z = 0 for r = 0 and r = d, together with the one-dimensional dis-

Fig. 11. The dimensionless pffiffiffiffiffiffi radial displacement Ur at the outer surface e for (from top to bottom) Dt=d e 2 ¼ 0:1, 1.0, versus the variable z= Dt 10.0, and 100.0. The solid curves indicate the numerical calculation of the general solution. The corresponding approximate solution in the e 2  1, is shown as a dashed curve which tends to thin limit, Dt=d the value
p/6 at z = 0; this curve is barely distinguishable from the envelope of the solutions.

pffiffiffiffiffiffi e =d. Thepone-dimenplacement, as a function of Dt ffiffiffiffiffiffi e =d. The sional displacement varies linearly with Dt axial displacements at r = 0 and r = d also show linear 5 behavior pffiffiffiffiffiffi for small times, but for t = 4.86 · 10 the value
2 e =d ¼ 3:7  10 of Dt is large enough that the displacements have begun to deviate from the onedimensional displacement, though the displacements at r = 0 and r = d are still approximately the same. In the limit of large times (not shown in Fig. 12), the difference between the axial displacements at r = 0 and r = d becomes small compared to their magnitude, and both approach one-third of the size of the one-dimensional displacement.

2006

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

Fig. 12. The dimensionless axial displacement Uz at z = 0 for r = 0 (solid curve) and r = d (dashed curve). The one-dimensional solution is also shown as the dotted line.

Fig. 13.ffi The dimensionless radial displacement Ur (solid curves) versus pffiffiffiffiffi e at r = d for (from top to bottom) Dt=d e 2 ¼ 0:001, 0.005, 0.01, z= Dt and 0.05. The corresponding dashed curves represent the approximae 2  1. tion obtained in the thick limit Dt=d

Fig. 13 shows the dimensionless radial pffiffiffiffiffidisplacement ffi e for a range Ur versus the scaled axial coordinate z= Dt e 2 . The corresponding of dimensionless times Dt=d e 2  1 are shown approximate solutions in the limit Dt=d 2 e as dashed curves. For Dt=d ¼ 0:001 the two curves are graphically indistinguishable, and the differences are just e 2 ¼ 0:005; for Dt=d e 2 ¼ 0:05 the noticeable for Dt=d deviation becomes significant. For the experimental time of t = 4.86 · 105 s, the dimensionless time of e 2 ¼ 1:4  10
3 is small enough that the radial disDt=d placement is still within the asymptotic regime.

4. Discussion The present work isolates the effects of different intrinsic diffusivities of the two species in a binary alloy and neglects the effect of differences in partial molar volume on deformation during diffusion. Clearly if the density

depends on composition, a sample will change shape as interdiffusion occurs. This effect can be added to the model and the diffusion equation will contain a term proportional to the product of the difference in partial molar volumes of the species and the gradient of the trace of the stress tensor. Thus the displacement field cannot be computed directly from the diffusion field independent of the moduli and the concentration field can not be computed a priori as was done in the present paper. Relaxation of the incompressibility condition would also require a different solution method that would require attention to the trace of the stress tensor. Intrinsic diffusion coefficients are often determined experimentally by measuring solute distribution in a diffusion couple, and measuring the Kirkendall plane location with respect to the Matano plane (defined by an equal area construction of the solute distribution [17]). As long as measurements are made on samples where e 2  1, and the concentration and marker shift is Dt=d measured in the center section of the sample, this procedure incurs only small errors that could be addressed through the results of this paper. We have computed a 3-D displacement field for a 1-D concentration profile and obtained a reasonable match with experimental data for the external shape of a diffusion couple. In the spirit of a small strain theory, we have ignored the axial and lateral distortion of the concentration profile due to the deformation. We find from the calculations that the axial strain reaches 0.05, suggesting that it would be useful to compare the present results with a finite strain analysis. In such an analysis, the diffusion equation (and the stress equilibrium equation) would have different forms in the reference configuration and the actual configuration. In this context the present paper solves the diffusion equation and applies the zero traction condition in the reference configuration in order to be able to use the 1-D error function solution and make the deformation predictions presented here. Cahn and Larche [18] and Mullins and Sekerka [19] also treat diffusion and deformation, but for the case of a network solid (a perfect lattice). In network models the strain is defined as a deformation of the lattice. In the Cahn and Larche description, the production and annihilation of vacancies occurs only at imperfections (grain boundaries, external surfaces, dislocations) separating perfect network regions. Boundary conditions on each of the individual perfect regions would then be used to respond to the requirements for vacancy production or annihilation (such an approach is also proposed to model creep). Although perfectly physical, numerical implementation of such an approach would require intense computations to treat more than a few defects. To our knowledge this sort of computation has not been performed. The present approach smears out this detail by assuming that the distance between defects is small compared to the diffusion distance and that any point

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

in the continuum can act as a source or sink of vacancies. Strain in this model is defined in terms of deformation of such a coarse-grained imperfect lattice, as might be described by the relative motion of inert markers embedded in the sample. We note that Adda and Philibert [20] consider balance equations that treat diffusion on networks where the number of sites is not conserved. They however do not treat deformation explicitly. In our analysis, we have assumed that D0 and the e are constant. From Fig. 3 we interdiffusion constant D see that as the atomic fraction of silver varies from XB = 0 to XB = 0.25, the interdiffusion constant varies e ¼ 5:7  10
14 m2/s to D e ¼ 7:9  10
14 m2/s, from D and D0 varies from D0 =
2.1 · 10
14 m2/s to D0 =
4.2 · 10
14 m2/s. Some of the discrepancy between the observed and calculated displacement fields shown in Fig. 6 may be due to this approximation. For example, the observed displacement field decays much more rapidly with z than the computed solution, and this may be due in part to using the diffusivities for XB = 0.125 to compute the far-field displacement field rather than using the diffusivities for XB = 0 and XB = 0.25. p Since ffiffiffiffiffiffi the solute field far from the interface e , the rate of decay of the displacement decays as Dt field must be dependent to some extent on the far-field e though it is unlikely that this would account values of D, entirely for the difference between the calculated and observed decay rates. In principle, the present analysis could be extended to include the effects of concentration-dependent diffusivities. A self-similar concentration profile may still be computed in this case, so that the velocity profiles could be obtained by Fourier analysis using the computed solute profiles and diffusivities. In theppresent work,p the ffiffiffiffiffiffi ffiffiffiffiffiffi computed bulge location e for small Dt e =d (see Fig. 13), and tends scales as Dt pffiffiffiffiffi ffi e =d (see Fig. 10). In the to a finite limit for large Dt experiment of Voigt and Ruth, the sample is much thicker than the diffusion p distance, and the bulge locaffi tion is found to scale with t [6]. In the limit of a sample much thicker than the diffusion distance, our model recovers the classical unidirectional Kirkendall displacement along the central axis for both slabs and rods. Near the surface of the sample, however, there are both axial and transverse displacements that do not exhibit one-dimensional behavior. The transverse displacement associated withp the ffiffiffiffiffiffisurface bulge penetrates to a depth e , and the axial position of the bulge of the order p offfiffiffiffiffiffi Dt e scales with Dt as well. Similarly, the axial displacement shows a variation near the surface with the same length scale. The occurrence pffiffiffiffiffiof ffi inhomogeneities near e with distance into the the surface that decay as Dt sample can be viewed in some sense as an example of St. Venant
s principle [21]; for thick samples there is no other length scale in the problem since the moduli do not affect the displacement field.

2007

5. Conclusions 1. A small-strain model employing a isotropic stress-free strain rate proportional to the difference in intrinsic (lattice) diffusivities, an incompressible elastic/plastic constitutive law, and a traction-free sample provides a good description to the external shape change data of Ruth et al. 2. The model also describes the full (internal) displacement field produced by one dimensional diffusion. Limiting analytical results are given for thin and thick pffiffiffiffiffi ffi e ). samples (compared to the diffusion distance Dt 3. For thick samples, the one dimensional Kirkendall displacement is recovered along the central axis for both slabs and rods. The axial displacement differs significantlypfrom ffiffiffiffiffiffi the one-dimensional result within e of the exterior surface. a distance Dt 4. For thin samples the stress vanishes and the displacement in the axial direction is one-third of the 1-D result for rods and one-half of the 1-D result for slabs. 5. A more comprehensive treatment could eliminate some of the simplifying assumptions in the model, allow an assessment of cumulative errors resulting from our simplified treatment, and allow a broader test for other alloy systems. An improved model might include the treatment of different partial molar volumes, small dislocation densities, concentrationdependent diffusion constants, void formation, multiphase systems, large deformation, more complex plasticity treatments, or relaxing the assumption of a one-dimensional diffusion field.

Acknowledgements The authors are grateful for helpful conversations with T.J. Burns, J.W. Cahn, C.A. Campbell, J.A. Dantzig, and D. Josell. C.A. Campbell kindly provided the diffusion analysis using DICTRA shown in Fig. 3. GBM, SRC and RFS are grateful for support of the Microgravity Research Division of NASA.

References [1] Philibert J. Atom movements: diffusion and mass transport in solids, translated from the French by Steven J. Rothman, Les Ulis, France: Editions de Physique, c1991. Monographies de physique. [2] da Correa Silva LC, Mehl RF. Interface and marker movements in diffusion in solid solutions of metals. Trans AIME 1951;191:155. [3] Huntington HB, Grone AR. Current-induced marker motion in gold wires. J Phys Chem Solids 1961;20:76. [4] Voigt R, Ruth V. A continuum theory of non-diffusional mass flow and arising stress fields in unbalanced diffusion systems. Defect Diffusion Forum 1997;143–147:483–8.

2008

W.J. Boettinger et al. / Acta Materialia 53 (2005) 1995–2008

[5] Busch P, Ruth V. The phenomenon of lateral mass flow associated with the Kirkendall effect at the gold-silver system. Z Metallkde 1989;80:238. [6] Voigt R, Ruth V. The development of surface contour changes during interdiffusion in silver gold alloy diffusion couples. J Phys – Condensed Matter 1995;7:2655. [7] Ruth V, Voigt R. Interdiffusion related internal volume changes and their effects on Kirkendall shift and surface contours of diffusion couples. Scripta Mater 1998;39:631. [8] Stephenson GB. Deformation during interdiffusion. Acta Metall Mater 1988;36:2663. [9] Happel J, Brenner H. Low Reynolds number hydrodynamics. Englewood Cliffs, NJ: Prentice-Hall; 1965. [10] Erdelyi A, Magnus W, Oberhettinger F, Tricomi FG. Tables of integral transforms, vol. I. New York: McGraw Hill; 1954. [11] SSOL2 Solutions Database (Version 2.1, 1999/2002), Scientific Group Thermodata Europe. [12] MOB2 Alloy Mobility Database (Version 2.0, 1999), ThermoCalc Software, Stockholm, Sweden. [13] Andersson JO, Helander T, Ho¨glund LH, Shi PF, Sundman B. THERMO-CALC & DICTRA, computational tools for materials science, Calphad 26 273, 2002.

[14] Abramowitz M, Stegun IA. Handbook of mathematical functions. New York: Dover; 1972. [15] Landau LD, Lifshitz EM. Theory of elasticity. New York: Pergamon Press; 1970. p. 75. ¨ berhuber CW, Kahaner [16] Piessens R, de Donker-Kapenga E, U DH. QUADPACK a subroutine package for automatic integration. Springer Series in Computational Mathematics, vol. 1. New York: Springer-Verlag; 1983. Subroutine QAWF is available through the NETLIB Repository atAvailable from: http:// www.netlib.org/ . [17] Glicksman ME. Diffusion in solids. New York: John Wiley & Sons; 2000. p. 306. [18] Larche FC, Cahn JW. The effect of self-stress on diffusion in solids. Acta Metall 1982;30:1835. [19] Mullins WW, Sekerka RF. On the thermodynamics of crystalline solids. J Chem Phys 1985;82:5192. [20] Adda Y, Philibert J. La Diffusion dans les solides, Saclay, Institut national des sciences et techniques nucleaires. Paris: Presses universitaires de France; 1966. [21] Fung YC. Foundations of solid mechanics. New Jersey: Prentice-Hall; 1965. p. 300.