- Email: [email protected]
Vacancy-driven stress relaxation in layers
Available online at www.sciencedirect.com
Acta Materialia 57 (2009) 4649–4657 www.elsevier.com/locate/actamat
Vacancy-driven stress relaxation in layers J. Svoboda a, F.D. Fischer b,* a
Institute of Physics of Materials, Academy of Sciences of the Czech Republic, Zˇizˇkova 22, CZ-616 62 Brno, Czech Republic b Institute of Mechanics, Montanuniversita¨t Leoben, Franz-Josef-Strasse 18, A-8700 Leoben, Austria Received 2 April 2009; received in revised form 5 June 2009; accepted 14 June 2009 Available online 21 July 2009
Abstract The role of vacancies together with their generally non-ideal sources and sinks has been investigated previously in the context of multicomponent diffusion. This concept is applied to a thin layer on a substrate subjected to an initial eigenstress state. The surface of the layer is assumed to act as an ideal source and sink for vacancies. The interface may act either as an ideal or as no source and sink for vacancies. Non-ideal sources and sinks are supposed in the bulk. The formation of a vacancy site fraction profile across the layer due to active sources and sinks and due to diffusion is studied, provoking a relaxation of the residual eigenstress state. A set of two nonlinear partial differential equations for both the vacancy site fraction and the eigenstress distributions is developed. The problem is reformulated in dimension-free quantities. The relation between the diffusivity and the activity of sources and sinks for vacancies in the bulk is characterized by a vacancy intensity parameter k. Numerical solutions are presented, and the influence of nonlinear terms in the evolution equations is evaluated. The simulations provide the evolution of the vacancy site fraction and of the hydrostatic stress profiles in the layer. The development of a film of newly deposited or removed atoms at the layer surface and of the total thickness of the layer is also calculated. The results of simulations are discussed. Ó 2009 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Substitutional diffusion; Vacancies; Annihilation; Residual stresses; Non-equilibrium processes
1. Introduction The compensation or relaxation of a stress state in a layer can be performed by several mechanisms (see e.g. Ref. [1] as a source of information in the literature). However, to our knowledge, a vacancy-driven stress relaxation in a layer is still an open problem. Significant contributions exist for very thin layers, consisting of only one grain in the thickness direction, with vacancy-driven diffusion along the grain boundaries (see Gibbs’ work from 1966 [2], and the works by Gao and co-workers [3]). A further important, yet unexplored, aspect is to study layers, consisting of grains with sources and sinks for vacancies (jogs at climbing dislocations) inside the grains. In addition, the surface of the layer and its usually incoherent interface to the substrate may act as ideal source and sink for vacancies. For *
Corresponding author. E-mail address: [email protected] (F.D. Fischer).
the case of a coherent interface one can assume no sources and sinks for vacancies existing there. The goal of this paper is to demonstrate the diffusive interaction of a single phase layer, as it is sketched in Fig. 1, with non-ideal sources and sinks for vacancies in its bulk with the surface acting as ideal source and sink for vacancies and with the interface acting either as ideal or no source and sink for vacancies. The relaxation of the initial eigenstress state and the development of the vacancy site fraction profile in the layer are investigated. 2. Problem formulation 2.1. Definitions, kinematics and conservation laws Let us assume that the material of the layer is characterized by a single component A and by vacancies, denominated by V. First of all it should be stated that an atomistic view of vacancies (see e.g. Ref. [4]) and a
1359-6454/$36.00 Ó 2009 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2009.06.016
4650
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
ð5Þ e_ gc ¼ a; a : d ¼ a; d is the unity tensor: This rather simple relation for e_ gc considers only the role of vacancy generation/annihilation and has no contribution due to volumetric misfit of the components A and V, since XA ¼ XV (for details see Ref. [5]). Finally it should be mentioned that treatments of nonideal sources and sinks hardly can be found. Here, the work by Garikipati et al. [8] should be given special mention, where an expression ðy V y eq V Þ=sG with sG being a time parameter is used for a. However, also ‘‘mixed” concepts are published as Ref. [9], using the case a = 0 and y V ¼ y eq V within one and the same concept. 2.2. Chemical potentials
Fig. 1. Geometrical setting and some related quantities.
mechanistic (or phenomenological) view of the vacancies exist. We follow the mechanistic view of vacancies enabling the diffusive motion of substitutional elements as e.g. outlined in Ref. [5]. Although the literature concerning multi-component diffusion is nearly numberless, only few treatments exist dealing in a rigorous way with vacancies. A main feature of vacancies is based on the fact that it is a non-conserving component and so it is advantageous to work with site fractions, in our case yA and yV, yielding ð1Þ y A þ y V ¼ 1; y V y A : Following the treatment by Svoboda et al. [5], conservation relations exist as y_ A ¼ Xdiv j A ay A ;
y_ V ¼ Xdiv j V þ ay A ;
ð2Þ
satisfying the balance relation ð3Þ j A þ j V ¼ 0: The rates are marked by dots, the vectors jA and jV are the diffusive fluxes of A and V, respectively, the quantity a is the rate, at which vacancies are generated (a > 0) or annihilated (a < 0). The volume corresponding to one mole of lattice positions is X. The same partial molar volumes of A, XA, and of V, XV , yield X ¼ XA ¼ XV . Two limits for a can be reported: No sources and/or sinks for vacancies: a = 0. Ideal sources and/or sinks for vacancies: a well-established feature exists addressing yV to its equilibrium eq value y eq V with the rate y_ V ¼ 0 under isothermal conditions. Then a can be found from Eq. (2) as a ¼ aeq ¼
X XA div j V ¼ div j A : yA yA
ð4Þ
This case is denoted as ‘‘Darken-case”, since Darken tacitly used ideal sources and sinks in his treatment of the Kirkendall problem (see e.g. Refs. [6,7]). However, the reader will recognize from Section 4.5 of this paper that Eq. (4) is only valid in a stress-free body. The vacancy generation/annihilation provokes a change in the deformation state. Consequently, we have the following eigenstrain rate e_ gc :
Let lA be the chemical potential of A and lV be the chemical potential of V, assumed to be lV ¼ Rg T ln ðy V =y eq V Þ with Rg being the gas constant and T the temperature supposed to be fixed. Finally the chemical part of the Gibbs energy per mole of lattice sites follows as ð6Þ gchem ¼ y A lA þ ð1 y A ÞlV : Consequently, one finds generalized chemical potentials, which can be formulated in the notation of this case (see Ref. [5], Eqs. (31), (32)) as @u lA ¼ lA lV þ umech ðXA XV Þ þ mech X @y A rH ðXA XV Þ ð7Þ @u ð8Þ lV ¼ d lV þ umech XV mech y A X rXV : @y A Note that lA is a scalar and lV a tensorial entity. The quantity umech represents the elastic strain energy (in joule per unit volume) and can be formulated as y A u0mech with u0mech as strain energy in a body with no vacancies (then @umech =@y A ¼ u0mech ); r is the stress tensor and rH the hydrostatic stress, d:r = 3rH. An equivalent formulation to Eqs. (7) and (8) was presented by Wu [10,11], and later by Varias and Massih [12], including the elastic strain energy term, and by Larche` and Voorhees [13] without strain energy terms. Recently a quite general formulation of the chemical potentials was presented in Ref. [14], but without any relations to the origin of the stress and strain terms in the potentials. A reformulation with the molar volume XA and u0mech yields lA ¼ lA lV þ u0mech XA ; ð9Þ ð10Þ lV ¼ dlV rXA ; lV ¼ lV : d=3 ¼ lV rH XA : Several papers have been published, mainly based on the ‘‘Darken-case”, where the term rHXA is added to Eq. (9) (see e.g. Ref. [15]), and Eq. (10) is ignored. 2.3. Constitutive equation for the fluxes Diffusion equations are derived by engaging Onsager’s maximum dissipation principle and Manning’s macroscopic theory of diffusion in Ref. [5] (Section 4). In the
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
4651
notation of this paper, the flux jA follows from Ref. [5] (Eqs. (40)–(46)) as
y_ V ¼ XA div j V þ a
y y Deq j A ¼ eqV A A grad lA : fy V XA Rg T
Deq A ½ð1 y V Þgrad y V fy eq V XA XA 0 ; ð17Þ u ð1 y V Þy V grad Rg T mech Deq divj V ¼ eqA fy V XA 2 3 D y V 1 RXgAT u0mech ðgrad y V Þ2 þ 6 7 5: ð18Þ 4 XA 0 XA 0 u Dy V þ grad y V grad Rg T umech Rg T mech
ð11Þ
Deq A is the tracer diffusion coefficient of A for the equilibrium concentration of vacancies, and f is the geometrical correlation factor (f = 0.7815 for face-centered cubic (fcc) and 0.7272 for body-centered cubic (bcc) alloys). Since lA, the chemical potential of phase A, can be considered as a constant quantity, or in other words, lA is (nearly) unaffected by the site fraction of vacancies with the consequence grad lA = 0, the vector grad lA in Eq. (9) can be calculated with lV ¼ Rg T lnðy V =y eq V Þ as grad lA ¼
Rg T grad y V þ XA gradðu0mech Þ: yV
ð12Þ
Using Eqs. (3), (11), and (12) yields j A ¼ j V ¼
y A Deq y V y A Deq A A gradðu0mech Þ: grad y V fy eq fy eq V XA V Rg T
ð13Þ
2.4. Constitutive equation for the vacancy-related eigenstrain rate a The simplest possible equation for the tensor a is, according to Eq. (54) in Ref. [5], a¼
1 l d; 3KXA V
ð14Þ
where K is the bulk viscosity factor, which can be related to the presence of jogs at dislocations being supposed to be the only sources and sinks for vacancies in the bulk of the layer. Furthermore, it is assumed that the Burgers vectors of dislocations are randomly oriented, and the jogs are homogeneously distributed in the layer. Thus the properties of the layer are assumed to be homogeneous and isotropic. The jogs at dislocations can be considered as ideal sources and sinks for vacancies. The Gibbs energy is dissipated by the transport of vacancies to or from the jogs by diffusion. The relation between the jog density and the value of the bulk viscosity factor K is given by Eq. (62) in Ref. [5]. Using Eq. (10) one has with Eqs. (5) and (14) 1 lV rH d ¼ e_ gc d; e_ gc ¼ a ¼ 3K XA 1 lV rH ¼ 3_egc : ð15Þ a¼ K XA 3. Evolution equations 3.1. General equation for the evolution of the vacancy distribution By assuming yA = 1 yV 1, one can rewrite Eq. (2) as
ð16Þ
and Eq. (13) as jV ¼
The symbol ‘‘D” represents the Laplace operator. Combining Eq. (18) with Eq. (15) and using lV ¼ Rg T ln ðy V =y eq V Þ yield for Eq. (16), 2 3 2 XA 0 eq D y 1 u Þ þ ðgrad y V V Rg T mech D 6 7 5 y_ V ¼ Aeq 4 XA XA fy V 0 0 u Dy V þ grad y V grad Rg T umech Rg T mech 1 Rg T yV : ð19Þ rH ln eq þ yV K XA Eq. (19) includes as unknown quantities yV and rH but also u0mech . Generally written, one has to know the full stress tensor, which enforces a coupled solution of the diffusion problem and the deformation (strain and stress) problem. In the next section we will specialize to a thin layer, which allows both an explicit calculation of the stress state and consequently also u0mech , so we can work finally with the two variables yV and rH. 3.2. Evolution equations of the vacancy distribution and eigenstress state in a thin layer Let us consider a thin layer of thickness d. The origin of the z-axis is positioned in the surface with the z-axis looking into the layer, 0 6 z 6 d. The vacancy site fraction yV as well as the hydrostatic stress rH are then functions of z and t; yV = yV(z, t), rH = rH(z, t). Consequently Eq. (19) can be simplified as 2 3 2 2 XA @ 0 eq y 1 u @y@zV þ 2 V D Rg T mech @z 5 y_ V ¼ Aeq 4 XA @ @y V 0 fy V u mech Rg T @z @z 1 Rg T yV þ : ð20Þ rH ln eq yV K XA The layer is assumed to be elastic, with E its Young’s modulus, v its Poisson’s ratio and to be fully constrained by the substrate. Furthermore, a small strain setting is applied with a linear decomposition of the strain tensor as e ¼ ee‘ þ egc . Since the components exx = eyy are fixed to a value e0, e.g. due to the manufacturing process and due e‘ to the constraint, one has ee‘ xx ¼ eyy ¼ e0 egc . The corresponding actual plane stress state components are
4652
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
Eegc ; rxx ¼ ryy ¼ r0 1v 2 Eegc : r0 ¼ 1v 3
rzz ¼ 0;
2 rH ¼ rxx 3
j
ð21Þ
Where r0 = 3/2rH0 is the initial in-plane homogeneous eigenstress state. The elastic strain component ezz follows from rzz = 0 as the sum of the elastic component ee‘ zz ¼ 2vegc =ð1 vÞ and egc as 1þv egc : ð22Þ 1v Since ezz is ow/oz with w being the according displacement in the z-direction, the displacement w(z, t) can be found with the boundary condition w(d) = 0 as Z Z d 1þv d ezz d~z ¼ egc ðz;~tÞd~z: ð23Þ wðz; tÞ ¼ 1v z z ezz ¼
The elastic strain energy can be expressed by rH as 9ð1 vÞ 2 u0mech ¼ rH : ð24Þ 4E Differentiation of Eq. (21), combined with Eq. (15), and using lV ¼ Rg T lnðy=y eq V Þ, yield 2E Rg T lnðy V =y eq V Þ ð25Þ rH : r_ H ¼ 9Kð1 vÞ XA Eqs. (20) and (25) represent a set of differential equations describing the system evolution with the state variables yV(z, t) and rH(z, t). 3.3. Boundary conditions at the surface and interface For the layer (0 6 z 6 d), the surface (z = 0) is assumed to be acting as ideal source and sink for vacancies. Furthermore, one can expect the interface (z = d) to be either an incoherent interface, acting as ideal source and sink for vacancies, or to be a coherent interface, acting as no source and sink for vacancies. The chemical potential of the vacancies is a tensorial entity lV in the layer, see Eq. (8). At the surface one can construct its projection lV m with m being the normal vector to the surface with components (0, 0, –1), see also Fig. 1. This projection can be interpreted as a thermodynamic traction vector, in analogy to a traction vector according to a stress tensor r. As we assume that the surface acts as an ideal source and sink for vacancies, the thermodynamic traction vector must be a zero vector. Thus, using Eq. (10) one has lV m ¼ lV r m ¼ 0:
ð26Þ
In the case of an interface acting as ideal source and sink for vacancies, the interface behaves like the surface, and the system is symmetric with respect to the plane z = d/2. Thus, due to the symmetry condition,
jV ;z ¼ 0
j
d at z ¼ : 2
ð28Þ
In the case of a coherent interface, acting as no source and sink for vacancies, the vacancy flux must disappear at the interface and thus
jV ;z ¼ 0
at z ¼ d:
ð29Þ
4. Results of simulations and their discussion 4.1. Solution of the differential equations The partial differential equations (20) and (25) can be transferred in a dimension-free form, introducing a dimension-free time parameter s: s¼
Deq A t 2 fy eq V d eff
ð30Þ
with t being the actual time. The quantity deff stands for deff = d in the case of an interface with no sources and sinks for vacancies and for deff = d/2 in the case of an interface with ideal sources and sinks for vacancies. Then one can introduce dimension-free coordinates (x/deff, y/deff, z/deff) ? (n, g, f). It is also advantageous to introduce the following dimension-free quantities and abbreviations: ~H ¼ r ~v ¼
rH ; E
~0 ¼ r
9ð1 vÞ ; 2
r0 ; E
j¼
~y ¼ y V =y eq V :
EXA ; Rg T
k¼
fRg Td 2eff ; KXA Deq A ð31Þ
Then the evolution equations (20) and (25) can be rewritten as ~2H d~y @ 2 j~vr ~y 1 ¼ 2 ds @f2 2 j~v @ @~y @~y ~2H y eq þ þ k ðj~ rH lnð~y ÞÞ; ð32Þ r V 2 @f @f @f d~ rH ky eq rH Þ: ¼ V ðln~y j~ ð33Þ ~v ds The boundary conditions follow from Eqs. (27)–(29), applied to Eq. (17) with 1 yV 1, as
ð27Þ
@~y j~v @ 2 ~y ð~ r Þ ¼ 0 for f ¼ 1: @f 2 @f H ð34Þ
For the interface we have two possible boundary conditions, expressed in jV,z being the z-component of the vector jV:
The set of the two partial differential equations (32) and (33), together with the boundary conditions (34), is of the so-called ‘‘spatially inhomogeneous reaction–diffusion
Assuming a stress-free surface, rxz = 0, ryz = 0, rzz = 0, relation (26) leads to lV = 0, yielding with lV ¼ Rg T lnðy V =y eq V Þ y V ¼ y eq V
at z ¼ 0:
~y ¼ 1 for f ¼ 0
and
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
(a)
1.010
1.008
3
2
4 5
1.006
yV /y V eq
type” (see e.g. Ref. [16]), and still a topic of ongoing research in applied mathematics; especially the long term asymptotic behaviour is of interest. The problem is, however, more complicated as one would assume at a first glance, since the associated stress field r, reflected by rH, is also not known in advance and evolves together with ~y due to the eigenstrain rate represented by e_ gc ¼ a, see Eq. (5). Since nonlinear terms appear in (32) and (33) as well in the boundary conditions (34), one cannot expect analytical solutions. Therefore, the equations (32) and (33), together with boundary conditions (34), are solved numerically by means of finite difference method on a grid of 101 equidistant points. The time integration is performed by Euler method.
4653
6 7 8 9
1.004
1.002
1.000
1
0.0
0.2
0.4
0.6
0.8
1.0
z/deff
4.2. Simulations of yV and rH profiles
(b)
The simulations are performed for the following values of parameters:
0.0010
1 2
0.0008
3
j
Values, typical for metals, are addressed to 9 v ¼ 3:15Þ. y eq V ¼ 10 ; j ¼ 10 and v ¼ 0:3ð~ The parameter k, which can be denominated as a vacancy intensity factor, relating the diffusion coefficient Deq A and the bulk viscosity factor K, is considered as a free parameter. In the simulations the values k = 1000, 10, 1 and 0.1 are chosen.
~H are chosen as The initial values of ~y and r ~y 0 ¼ 1 and r ~H 0 ¼ 103 , respectively. Four sets of diagrams, Figs. 2–5, are presented according to four distinct values of the vacancy intensity factor k, which obtains k ? 1 for ideal sources and sinks and k ? 0 for no sources and sinks in the bulk. Figs. 2–5 demonstrate the relative vacancy site fraction ~y ¼ y V =y eq V profile across the layer, starting from the surface towards z = deff. A common feature of all the four diagrams is that ~y changes quickly from its initial state (curve label 1), ~y ¼ 1, to its highest values (curve label 2) during the first time increment. Such a supersaturation with vacancies occurs due to the generation of vacancies at sources in the bulk induced by the stress state. This process is studied in detail with much smaller time increments in Fig. 6 for k = 10. One observes that after approximately 20 small time increments, which means approximately for s ¼ 0:4 (a time period much shorter than one distinct time increment in Figs. 2–5) an envelope of the profiles is formed, pointing to a maximum value of ~y ¼ 1:008. Obviously the bulk provides enough sources of vacancies, and the initial equilibrium distribution ~y ¼ 1 is quickly forgotten. However, with an ongoing diffusion process the vacancies move towards the surface, provoke further generation of vacancies causing stress relaxation and thus reduce the supersaturation in the bulk (see Figs. 2–5). After a certain time interval, say 20 time increments, the equilibrium profile ~y ¼ 1 is approached again. This coupled vacancy generation/annihilation and diffusion process in the bulk and the surface acting as ideal sink/ source for vacancies allows for a rather gradual reduction
4 0.0006
5
σH /E
j
6 7
0.0004
8 9
0.0002
0.0000 0.0
0.2
0.4
0.6
0.8
1.0
z/d eff 9 ~H profiles for y eq Fig. 2. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3ð~v ¼ 3:15Þ and k ¼ 1000. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 2.5 107 with respect to the preceding one.
of the eigenstress state, represented by the profile of relative ~H for all four values of k (see Figs. 2–5). hydrostatic stress r The eigenstress state decreases most quickly near the surface. But with ongoing time (curve labels >10) the eigenstress state is significantly relaxed and finally disappears nearly totally across the whole layer thickness. 4.3. The role of nonlinearities in the differential equations One can see immediately from Figs. 2–5 that ~y is close to 1, which allows a replacement of ~y by 1 þ D~y (the symbol D should not be confused with the Laplace Operator D) and a linearization of lnð~y Þ D~y . Then one can rewrite the set of differential equations (32) and (33), ordered with respect to linear and nonlinear terms, as dD~y @ 2 D~y ~H Þ; ¼ kðD~y j~ rH Þ þ nD~y ðD~y ; r ds @f2 d~ rH y eq rH Þ: ¼ V kðD~y j~ ~v ds
ð35Þ ð36Þ
4654
(a)
1.008
2
3
yV /y V eq
1.006
1.003
3
2
4
4 5
5 6
1.002
6 7
1.004
8 9
yV /y V eq
(a)
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
7 8 9
1.001
1.002
1.000
1.000
1 0.0
0.2
0.4
0.6
0.8
1 0.0
1.0
0.2
0.4
0.6
(b)
1
0.0010
1
2
2
3
0.0008
3
0.0008
4 5
5
0.0006
6 7 8
0.0004
6
0.0006
7
σH /E
σH /E
4
9
8
9
0.0004
0.0002
0.0002
0.0000
0.0000
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
9 ~H profiles for y eq Fig. 3. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3ð~v ¼ 3:15Þ and k ¼ 10. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 2.5 107 with respect to the preceding one.
The nonlinear contribution nD~y to the first differential equation (35) follows after rearrangement as ( ) 2 y j~v @D~y @ð~ r2H Þ @ 2 ð~ r2H Þ eq @D~ ~H Þ ¼ y V : þ þ nD~y ðD~y ; r @f 2 @f @f @f2 ð37Þ The corresponding boundary conditions (34) can be reformulated as D~y ¼ 0 for f ¼ 0 and
@D~y ¼ nr~H ¼ 0 for f ¼ 1 @f
0.6
0.8
1.0
z/deff
z/d eff
ð38Þ
with nr~H ¼ j~v
1.0
z/deff
z/deff
(b) 0.0010
0.8
9 ~H profiles for y eq Fig. 4. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3ð~v ¼ 3:15Þ and k ¼ 1. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 5.0 107with respect to the preceding one.
the boundary conditions (39) with simulation results without the nonlinear terms. The results of simulations differ within 0.2% to the simulations including all terms and cannot be distinguished in Figs. 2–6. If the nonlinear terms can be neglected, the differential equation (35) can be reformulated as d @ 2 D~y V~ ~H ¼ D~y þ eq r ; ds yV @f2 d~ rH k eq ¼ y V ðD~y j~ rH Þ ~v ds
ð40Þ ð41Þ
with the boundary conditions @ð~ r2H Þ : @f
ð39Þ
The role of the nonlinear contributions is checked numerically by comparing simulation results including the terms nD~y in the differential equation (35) and nr~H in
D~y ¼ 0 for f ¼ 0
and
@D~y ¼ 0 for f ¼ 1: @f
ð42Þ
Two coupled partial differential equations for the two vari~H ðf; sÞ are to be solved. The following ables D~y ðf; sÞ and r consequences can be drawn:
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
(a)
1.0004 j
2
1.0003
3
4
yV /y V eq
5 6 7
1.0002
8 9 j
1.0001
1.0000
1 0.0
0.2
0.4
0.6
0.8
1.0
z/d eff
(b)
0.0010
1
0.0008
~v 1 @ 2 D~y ¼ eq j~ rH ; yV k @f2 d~ rH y eq y eq rH ¼ k V D~y : þ k V j~ ~v ~v ds
0.0006
ð43Þ ð44Þ
This formulation allows in principle a semianalytic solution of Eq. (43) for D~y via the Green’s function to the operator on the left side and finally the construction of an integro~H . differential equation for r
3 4
σH /E
~H are replaced by D~y and r ~H in the set If D~y and r of differential equations (40) and (41) as well as in the boundary conditions (42), all equations remain unchanged. This means that a solution for the ‘‘tension case” ð~ rH 0 > 0Þ provides automatically the solution for the ‘‘compression case” ð~ rH0 < 0Þ. The quantity ð~v=y eq rH is usually strongly dominant in V Þ~ 9 comparison to D~y , because 1=y eq V is in order of 10 . Since ~ H are smooth, see Figs. 2–5, one the functions D~y and r can neglect D~y on the left side of Eq. (40) yielding by combination of Eqs. (40) and (41)
D~y
2
5 6 7
0.0004
8 9
0.0002
0.0000 0.0
0.2
0.4
0.6
0.8
1.0
z/deff 9 ~H profiles for y eq Fig. 5. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3 ð~v ¼ 3:15Þ and k ¼ 0:1. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 5.0 108 with respect to the preceding one.
7 6 5
1.006
4 1.004
1.002
0.0
0.2
0.4
0.6
0.8
Let h = h(t) be the temporal change of deff from its initial value and ~h ¼ h=d eff . Two contributions to ~ h can be observed, the first one due to the deposition of atoms at the surface, denoted as ~hv , and the second one only due to the deformation process, described by egc and the corresponding elastic strain state, denoted as ~hd . The rate h_ v follows from the mass balance at the surface, expressed in the reference configuration, see e.g. Ref. [17], Section 4.1, as y ð45Þ h_ v A ¼ jA;Z ð0; tÞ ¼ jV ;Z ð0; tÞ: XA
3
h_ v ¼ XA jV ;z ð0; tÞ
2
If a can be approximated by XA div jv, see Eq. (16), which is acceptable for times s > 0.4, (see Fig. 6), then integration of Eq. (16) and using Eq. (15) provide Z d eff Z d eff @jV ;z dz ¼ XA jV ;z ð0; tÞ ¼ 3 XA ð48Þ e_ gc dz: @z 0 0
1
1.000
4.4. Simulations of the evolution of the layer thickness
Inserting Eq. (17), together with Eqs. (30) and (31), into Eq. (45) provides d~hv y @ eq @~ 2 ~H Þ : ¼ y V ~y ðj~vr ð46Þ ds @f @f ~H in A direct relation between ~hv and h~ rH i, the average of r the interval 0 6 f 6 1, can be derived. Using yA 1 in Eq. (45) one obtains
1.008
yV /y V eq
4655
1.0
z/d eff 9 Fig. 6. Evolution of ~y profile for y eq V ¼ 10 ; j ¼ 10; v ¼ 0:3 ð~v ¼ 3:15Þ and k ¼ 10. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . ., 21 correspond to values of s increased by Ds = 0.02 with respect to the preceding one.
Then putting Eqs. (47) and (48) together gives Z d eff h_ V ¼ 3 e_ gc dz: 0
ð47Þ
ð49Þ
4656
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
Differentiation of Eq. (21)3 with respect to time provides 3ð1 vÞr_ H ; 2E
k=0.1
ð50Þ
and its insertion into Eq. (49) yields Z 9ð1 vÞ d eff 9ð1 vÞd eff _ hri: h_ V ¼ r_ H dz ¼ 2E 2E 0
0.0008
0.0006
ð51Þ
Integration of Eq. (51) in time and expression of the variables in dimension-free ones yield
<σH /E>
e_ gc
0.0010
~ ~H 0 Þ: hV ¼ ~vðh~ rH i r ð52Þ ~ The contribution hd due to the deformation process can directly be related to the displacement at the surface, see Eq. (23), as wð0; tÞ ~ : hd ¼ d eff
ð53Þ
Note the negative sign on the right side of Eq. (53), since w is counted as positive in the z-direction. Insertion of egc from Eq. (21) in Eq. (23) yields 3ð1 þ vÞ ~ hd ðtÞ ¼ ð~ rH 0 h~ rH iÞ: ð54Þ 2 It is interesting to note that both contributions to the hd are proportional change of the layer thickness, ~ hV and ~ 0 rH iÞ. Therefore, the reader is provided with a to ð~ rH h~ diagram for h~ rH i for various values of k in Fig. 7. The contributions of both deformation mechanisms are in a fixed hV ¼ ð1 þ vÞ=ð3ð1 vÞÞ. relation, depending only v as ~ hd = ~ 4.5. Limit cases for ideal and no sources and sinks for vacancies in the bulk For the layer with ideal sources and sinks for vacancies in the bulk the vacancy intensity parameter k tends to 1, rH must be limited for any s > 0. Then and the quantity d~ds it follows from Eq. (33) that ln ~y ¼ j~ rH , which is equivalent to lV ¼ 0, see Eq. (10). Its differentiation with respect to s rH 1 d~y provides j~ ¼ d~ds , and Eq. (33) can be rewritten as y ds ~v d~y ¼ kðln ~y j~ rH Þ: j~y y eq V ds
ð55Þ
Insertion of Eq. (55) with ln ~y ¼ j~ rH into Eq. (32) provides ~v ~vln2 ~y d~y @2 ~ ¼ y 1 1þ 2j j~y y eq ds @f2 V 2 ~v @ y y eq @~ 2 @~ þ : ð56Þ ln ~y yV 2j @f @f @f ~v @ @~y ~y ðln2 ~y Þ ¼ 0 for f ¼ d eff : 2j @f @f
Thus, the problem reduces to the only one partial differen~H is tial equation (56) with boundary conditions (57), and r
0.0002
k=10
0.0000
k=1000 2x108
4x108
τ
6x10 8
8x108
9
1x10
Fig. 7. Evolution of the relative hydrostatic stress h~ rH i as function of s for 9 v ¼ 3:15Þ and k ¼ 1000; 10; 1 the parameters y eq V ¼ 10 ; j ¼ 10; v ¼ 0:3 ð~ and 0:1. Note that the curve for k = 1000 de facto coincides for h~ rH i due to ideal sources and sinks.
rH i de facto coincides given by ln ~y ¼ j~ rH , whose average h~ with the case k = 1000 in Fig. 7. For no sources and sinks for vacancies in the bulk of the layer, the vacancy intensity parameter is 0. Since vacancies are neither generated nor annihilated, no change in the volume is possible (see the eigenstrain rate e_ gc ¼ 0 for a ¼ 0). ~H Eq. (33) yields then drH/ds = 0, which means that r remains constant in time. Only the vacancy site fraction profile evolves according to 2 ~2H d~y @2 j~vr j~v @ y y eq @~ 2 @~ ~H þ : ¼ 2 ~y 1 r yV 2 ds @f 2 @f @f @f ð58Þ The boundary conditions remain the same as given by Eq. (34). Again, the problem reduces to the solution of only one partial differential equation (58). 5. Conclusions The results achieved in this paper can be summarized as follows: j
j
j
ð57Þ
0.0004
0
The boundary conditions (34) can be then modified as ~y ¼ 1 for f ¼ 0 and
k=1
A model for simulation of stress relaxation in a layer on a substrate by means of diffusion and generation/annihilation of vacancies has been developed. The model assumes ideal sources and sinks for vacancies on the layer surface, ideal or no sources and sinks for vacancies at the interface and non-ideal sources and sinks for vacancies in the bulk of the layer. The driving force for diffusion and generation/annihilation of vacancies stems from the generalized chemical potentials taking into account the actual stress state according to Ref. [5]. The model leads to a boundary value problem, given by a set of two nonlinear partial differential equations for the vacancy site fraction and the hydrostatic stress. The equations are expressed in dimension-free quantities and solved numerically by finite differences.
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657 j
The results of simulations are presented for various relations between the diffusivity and activity of non-ideal sources and sinks for vacancies in the bulk. The analysis is completed by calculations of the change of the layer thickness due to the deformation state and due to the deposition/removing of a film of atoms at the layer surface.
6. Outlook The stress relaxation in surface layers can be performed by several mechanisms [1], and their coupling with the vacancy mechanism seems to be a challenge for the nearest future. One can imagine a coupling of the vacancy mechanism with diffusional and dislocation creep, which could complete the picture on stress relaxation in layers at elevated temperatures. Acknowledgements Financial support by Materials Center Leoben Forschung GmbH (project SP19) and by Research Plan of Institute of Physics of Materials (project CEZ:AV0Z20410507) is gratefully acknowledged. The cooperation of authors is also supported in the frame of COST action P19.
4657
References [1] Freund LB, Suresh S. Thin film materials, stress, defect formation and surface evolution. Cambridge: Cambridge University Press; 2003. [2] Gibbs GB. Philos Mag 1966;13:589. [3] Bhandakkar T, Wei Y, Gao H. Math Mech Solids 2009;14:179. [4] Lindgren LE, Domkin K, Hansson S. Mech Mater 2008;40:907. [5] Svoboda J, Fischer FD, Fratzl P. Acta Mater 2006;54:3043. [6] Svoboda J, Fischer FD, Gamsja¨ger E. Acta Mater 2008;56:351. [7] Suo Z, Kubair DV, Evans AG, Clarke DR, Tolpygo VK. Acta Mater 2003;51:959. [8] Garikipati K, Bassman L, Deal M. J Mech Phys Solids 2001;49:1209. [9] Tsurkov I, Grich WM, Gross TS. Analysis of diffusional stress relaxation in submicron Cu interconnect structures using the model with enhanced vacancy diffusivity in grain boundary region. In: Brebbia CA, editor. High performance structures and materials III. Wessex: WIT Press; 2006. p. 85. [10] Wu CH. J Mech Phys Solids 2001;49:1771. [11] Wu CH. Int J Fract 2007;147:227. [12] Varias AG, Massih AR. J Mech Phys Solids 2002;50:1469. [13] Larche` FC, Voorhees PW. Defect Diffus Forum 1996;129–130:31. [14] Spicer JB, Dikmelik Y. J Appl Phys 2008;104:023528-1-6. [15] Danielewski M, Wierzba B, Bachorczyk-Nagy R, Pietrzyk M. J Phase Equilibr Diff 2006;27:691. [16] Di Francesco M, Fellner K, Markowich PA. Proc R Soc A 2008;464:3273. [17] Fischer FD, Waitz T, Vollath D, Simha NK. Prog Mater Sci 2008;53:481.
Acta Materialia 57 (2009) 4649–4657 www.elsevier.com/locate/actamat
Vacancy-driven stress relaxation in layers J. Svoboda a, F.D. Fischer b,* a
Institute of Physics of Materials, Academy of Sciences of the Czech Republic, Zˇizˇkova 22, CZ-616 62 Brno, Czech Republic b Institute of Mechanics, Montanuniversita¨t Leoben, Franz-Josef-Strasse 18, A-8700 Leoben, Austria Received 2 April 2009; received in revised form 5 June 2009; accepted 14 June 2009 Available online 21 July 2009
Abstract The role of vacancies together with their generally non-ideal sources and sinks has been investigated previously in the context of multicomponent diffusion. This concept is applied to a thin layer on a substrate subjected to an initial eigenstress state. The surface of the layer is assumed to act as an ideal source and sink for vacancies. The interface may act either as an ideal or as no source and sink for vacancies. Non-ideal sources and sinks are supposed in the bulk. The formation of a vacancy site fraction profile across the layer due to active sources and sinks and due to diffusion is studied, provoking a relaxation of the residual eigenstress state. A set of two nonlinear partial differential equations for both the vacancy site fraction and the eigenstress distributions is developed. The problem is reformulated in dimension-free quantities. The relation between the diffusivity and the activity of sources and sinks for vacancies in the bulk is characterized by a vacancy intensity parameter k. Numerical solutions are presented, and the influence of nonlinear terms in the evolution equations is evaluated. The simulations provide the evolution of the vacancy site fraction and of the hydrostatic stress profiles in the layer. The development of a film of newly deposited or removed atoms at the layer surface and of the total thickness of the layer is also calculated. The results of simulations are discussed. Ó 2009 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Substitutional diffusion; Vacancies; Annihilation; Residual stresses; Non-equilibrium processes
1. Introduction The compensation or relaxation of a stress state in a layer can be performed by several mechanisms (see e.g. Ref. [1] as a source of information in the literature). However, to our knowledge, a vacancy-driven stress relaxation in a layer is still an open problem. Significant contributions exist for very thin layers, consisting of only one grain in the thickness direction, with vacancy-driven diffusion along the grain boundaries (see Gibbs’ work from 1966 [2], and the works by Gao and co-workers [3]). A further important, yet unexplored, aspect is to study layers, consisting of grains with sources and sinks for vacancies (jogs at climbing dislocations) inside the grains. In addition, the surface of the layer and its usually incoherent interface to the substrate may act as ideal source and sink for vacancies. For *
Corresponding author. E-mail address: [email protected] (F.D. Fischer).
the case of a coherent interface one can assume no sources and sinks for vacancies existing there. The goal of this paper is to demonstrate the diffusive interaction of a single phase layer, as it is sketched in Fig. 1, with non-ideal sources and sinks for vacancies in its bulk with the surface acting as ideal source and sink for vacancies and with the interface acting either as ideal or no source and sink for vacancies. The relaxation of the initial eigenstress state and the development of the vacancy site fraction profile in the layer are investigated. 2. Problem formulation 2.1. Definitions, kinematics and conservation laws Let us assume that the material of the layer is characterized by a single component A and by vacancies, denominated by V. First of all it should be stated that an atomistic view of vacancies (see e.g. Ref. [4]) and a
1359-6454/$36.00 Ó 2009 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actamat.2009.06.016
4650
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
ð5Þ e_ gc ¼ a; a : d ¼ a; d is the unity tensor: This rather simple relation for e_ gc considers only the role of vacancy generation/annihilation and has no contribution due to volumetric misfit of the components A and V, since XA ¼ XV (for details see Ref. [5]). Finally it should be mentioned that treatments of nonideal sources and sinks hardly can be found. Here, the work by Garikipati et al. [8] should be given special mention, where an expression ðy V y eq V Þ=sG with sG being a time parameter is used for a. However, also ‘‘mixed” concepts are published as Ref. [9], using the case a = 0 and y V ¼ y eq V within one and the same concept. 2.2. Chemical potentials
Fig. 1. Geometrical setting and some related quantities.
mechanistic (or phenomenological) view of the vacancies exist. We follow the mechanistic view of vacancies enabling the diffusive motion of substitutional elements as e.g. outlined in Ref. [5]. Although the literature concerning multi-component diffusion is nearly numberless, only few treatments exist dealing in a rigorous way with vacancies. A main feature of vacancies is based on the fact that it is a non-conserving component and so it is advantageous to work with site fractions, in our case yA and yV, yielding ð1Þ y A þ y V ¼ 1; y V y A : Following the treatment by Svoboda et al. [5], conservation relations exist as y_ A ¼ Xdiv j A ay A ;
y_ V ¼ Xdiv j V þ ay A ;
ð2Þ
satisfying the balance relation ð3Þ j A þ j V ¼ 0: The rates are marked by dots, the vectors jA and jV are the diffusive fluxes of A and V, respectively, the quantity a is the rate, at which vacancies are generated (a > 0) or annihilated (a < 0). The volume corresponding to one mole of lattice positions is X. The same partial molar volumes of A, XA, and of V, XV , yield X ¼ XA ¼ XV . Two limits for a can be reported: No sources and/or sinks for vacancies: a = 0. Ideal sources and/or sinks for vacancies: a well-established feature exists addressing yV to its equilibrium eq value y eq V with the rate y_ V ¼ 0 under isothermal conditions. Then a can be found from Eq. (2) as a ¼ aeq ¼
X XA div j V ¼ div j A : yA yA
ð4Þ
This case is denoted as ‘‘Darken-case”, since Darken tacitly used ideal sources and sinks in his treatment of the Kirkendall problem (see e.g. Refs. [6,7]). However, the reader will recognize from Section 4.5 of this paper that Eq. (4) is only valid in a stress-free body. The vacancy generation/annihilation provokes a change in the deformation state. Consequently, we have the following eigenstrain rate e_ gc :
Let lA be the chemical potential of A and lV be the chemical potential of V, assumed to be lV ¼ Rg T ln ðy V =y eq V Þ with Rg being the gas constant and T the temperature supposed to be fixed. Finally the chemical part of the Gibbs energy per mole of lattice sites follows as ð6Þ gchem ¼ y A lA þ ð1 y A ÞlV : Consequently, one finds generalized chemical potentials, which can be formulated in the notation of this case (see Ref. [5], Eqs. (31), (32)) as @u lA ¼ lA lV þ umech ðXA XV Þ þ mech X @y A rH ðXA XV Þ ð7Þ @u ð8Þ lV ¼ d lV þ umech XV mech y A X rXV : @y A Note that lA is a scalar and lV a tensorial entity. The quantity umech represents the elastic strain energy (in joule per unit volume) and can be formulated as y A u0mech with u0mech as strain energy in a body with no vacancies (then @umech =@y A ¼ u0mech ); r is the stress tensor and rH the hydrostatic stress, d:r = 3rH. An equivalent formulation to Eqs. (7) and (8) was presented by Wu [10,11], and later by Varias and Massih [12], including the elastic strain energy term, and by Larche` and Voorhees [13] without strain energy terms. Recently a quite general formulation of the chemical potentials was presented in Ref. [14], but without any relations to the origin of the stress and strain terms in the potentials. A reformulation with the molar volume XA and u0mech yields lA ¼ lA lV þ u0mech XA ; ð9Þ ð10Þ lV ¼ dlV rXA ; lV ¼ lV : d=3 ¼ lV rH XA : Several papers have been published, mainly based on the ‘‘Darken-case”, where the term rHXA is added to Eq. (9) (see e.g. Ref. [15]), and Eq. (10) is ignored. 2.3. Constitutive equation for the fluxes Diffusion equations are derived by engaging Onsager’s maximum dissipation principle and Manning’s macroscopic theory of diffusion in Ref. [5] (Section 4). In the
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
4651
notation of this paper, the flux jA follows from Ref. [5] (Eqs. (40)–(46)) as
y_ V ¼ XA div j V þ a
y y Deq j A ¼ eqV A A grad lA : fy V XA Rg T
Deq A ½ð1 y V Þgrad y V fy eq V XA XA 0 ; ð17Þ u ð1 y V Þy V grad Rg T mech Deq divj V ¼ eqA fy V XA 2 3 D y V 1 RXgAT u0mech ðgrad y V Þ2 þ 6 7 5: ð18Þ 4 XA 0 XA 0 u Dy V þ grad y V grad Rg T umech Rg T mech
ð11Þ
Deq A is the tracer diffusion coefficient of A for the equilibrium concentration of vacancies, and f is the geometrical correlation factor (f = 0.7815 for face-centered cubic (fcc) and 0.7272 for body-centered cubic (bcc) alloys). Since lA, the chemical potential of phase A, can be considered as a constant quantity, or in other words, lA is (nearly) unaffected by the site fraction of vacancies with the consequence grad lA = 0, the vector grad lA in Eq. (9) can be calculated with lV ¼ Rg T lnðy V =y eq V Þ as grad lA ¼
Rg T grad y V þ XA gradðu0mech Þ: yV
ð12Þ
Using Eqs. (3), (11), and (12) yields j A ¼ j V ¼
y A Deq y V y A Deq A A gradðu0mech Þ: grad y V fy eq fy eq V XA V Rg T
ð13Þ
2.4. Constitutive equation for the vacancy-related eigenstrain rate a The simplest possible equation for the tensor a is, according to Eq. (54) in Ref. [5], a¼
1 l d; 3KXA V
ð14Þ
where K is the bulk viscosity factor, which can be related to the presence of jogs at dislocations being supposed to be the only sources and sinks for vacancies in the bulk of the layer. Furthermore, it is assumed that the Burgers vectors of dislocations are randomly oriented, and the jogs are homogeneously distributed in the layer. Thus the properties of the layer are assumed to be homogeneous and isotropic. The jogs at dislocations can be considered as ideal sources and sinks for vacancies. The Gibbs energy is dissipated by the transport of vacancies to or from the jogs by diffusion. The relation between the jog density and the value of the bulk viscosity factor K is given by Eq. (62) in Ref. [5]. Using Eq. (10) one has with Eqs. (5) and (14) 1 lV rH d ¼ e_ gc d; e_ gc ¼ a ¼ 3K XA 1 lV rH ¼ 3_egc : ð15Þ a¼ K XA 3. Evolution equations 3.1. General equation for the evolution of the vacancy distribution By assuming yA = 1 yV 1, one can rewrite Eq. (2) as
ð16Þ
and Eq. (13) as jV ¼
The symbol ‘‘D” represents the Laplace operator. Combining Eq. (18) with Eq. (15) and using lV ¼ Rg T ln ðy V =y eq V Þ yield for Eq. (16), 2 3 2 XA 0 eq D y 1 u Þ þ ðgrad y V V Rg T mech D 6 7 5 y_ V ¼ Aeq 4 XA XA fy V 0 0 u Dy V þ grad y V grad Rg T umech Rg T mech 1 Rg T yV : ð19Þ rH ln eq þ yV K XA Eq. (19) includes as unknown quantities yV and rH but also u0mech . Generally written, one has to know the full stress tensor, which enforces a coupled solution of the diffusion problem and the deformation (strain and stress) problem. In the next section we will specialize to a thin layer, which allows both an explicit calculation of the stress state and consequently also u0mech , so we can work finally with the two variables yV and rH. 3.2. Evolution equations of the vacancy distribution and eigenstress state in a thin layer Let us consider a thin layer of thickness d. The origin of the z-axis is positioned in the surface with the z-axis looking into the layer, 0 6 z 6 d. The vacancy site fraction yV as well as the hydrostatic stress rH are then functions of z and t; yV = yV(z, t), rH = rH(z, t). Consequently Eq. (19) can be simplified as 2 3 2 2 XA @ 0 eq y 1 u @y@zV þ 2 V D Rg T mech @z 5 y_ V ¼ Aeq 4 XA @ @y V 0 fy V u mech Rg T @z @z 1 Rg T yV þ : ð20Þ rH ln eq yV K XA The layer is assumed to be elastic, with E its Young’s modulus, v its Poisson’s ratio and to be fully constrained by the substrate. Furthermore, a small strain setting is applied with a linear decomposition of the strain tensor as e ¼ ee‘ þ egc . Since the components exx = eyy are fixed to a value e0, e.g. due to the manufacturing process and due e‘ to the constraint, one has ee‘ xx ¼ eyy ¼ e0 egc . The corresponding actual plane stress state components are
4652
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
Eegc ; rxx ¼ ryy ¼ r0 1v 2 Eegc : r0 ¼ 1v 3
rzz ¼ 0;
2 rH ¼ rxx 3
j
ð21Þ
Where r0 = 3/2rH0 is the initial in-plane homogeneous eigenstress state. The elastic strain component ezz follows from rzz = 0 as the sum of the elastic component ee‘ zz ¼ 2vegc =ð1 vÞ and egc as 1þv egc : ð22Þ 1v Since ezz is ow/oz with w being the according displacement in the z-direction, the displacement w(z, t) can be found with the boundary condition w(d) = 0 as Z Z d 1þv d ezz d~z ¼ egc ðz;~tÞd~z: ð23Þ wðz; tÞ ¼ 1v z z ezz ¼
The elastic strain energy can be expressed by rH as 9ð1 vÞ 2 u0mech ¼ rH : ð24Þ 4E Differentiation of Eq. (21), combined with Eq. (15), and using lV ¼ Rg T lnðy=y eq V Þ, yield 2E Rg T lnðy V =y eq V Þ ð25Þ rH : r_ H ¼ 9Kð1 vÞ XA Eqs. (20) and (25) represent a set of differential equations describing the system evolution with the state variables yV(z, t) and rH(z, t). 3.3. Boundary conditions at the surface and interface For the layer (0 6 z 6 d), the surface (z = 0) is assumed to be acting as ideal source and sink for vacancies. Furthermore, one can expect the interface (z = d) to be either an incoherent interface, acting as ideal source and sink for vacancies, or to be a coherent interface, acting as no source and sink for vacancies. The chemical potential of the vacancies is a tensorial entity lV in the layer, see Eq. (8). At the surface one can construct its projection lV m with m being the normal vector to the surface with components (0, 0, –1), see also Fig. 1. This projection can be interpreted as a thermodynamic traction vector, in analogy to a traction vector according to a stress tensor r. As we assume that the surface acts as an ideal source and sink for vacancies, the thermodynamic traction vector must be a zero vector. Thus, using Eq. (10) one has lV m ¼ lV r m ¼ 0:
ð26Þ
In the case of an interface acting as ideal source and sink for vacancies, the interface behaves like the surface, and the system is symmetric with respect to the plane z = d/2. Thus, due to the symmetry condition,
jV ;z ¼ 0
j
d at z ¼ : 2
ð28Þ
In the case of a coherent interface, acting as no source and sink for vacancies, the vacancy flux must disappear at the interface and thus
jV ;z ¼ 0
at z ¼ d:
ð29Þ
4. Results of simulations and their discussion 4.1. Solution of the differential equations The partial differential equations (20) and (25) can be transferred in a dimension-free form, introducing a dimension-free time parameter s: s¼
Deq A t 2 fy eq V d eff
ð30Þ
with t being the actual time. The quantity deff stands for deff = d in the case of an interface with no sources and sinks for vacancies and for deff = d/2 in the case of an interface with ideal sources and sinks for vacancies. Then one can introduce dimension-free coordinates (x/deff, y/deff, z/deff) ? (n, g, f). It is also advantageous to introduce the following dimension-free quantities and abbreviations: ~H ¼ r ~v ¼
rH ; E
~0 ¼ r
9ð1 vÞ ; 2
r0 ; E
j¼
~y ¼ y V =y eq V :
EXA ; Rg T
k¼
fRg Td 2eff ; KXA Deq A ð31Þ
Then the evolution equations (20) and (25) can be rewritten as ~2H d~y @ 2 j~vr ~y 1 ¼ 2 ds @f2 2 j~v @ @~y @~y ~2H y eq þ þ k ðj~ rH lnð~y ÞÞ; ð32Þ r V 2 @f @f @f d~ rH ky eq rH Þ: ¼ V ðln~y j~ ð33Þ ~v ds The boundary conditions follow from Eqs. (27)–(29), applied to Eq. (17) with 1 yV 1, as
ð27Þ
@~y j~v @ 2 ~y ð~ r Þ ¼ 0 for f ¼ 1: @f 2 @f H ð34Þ
For the interface we have two possible boundary conditions, expressed in jV,z being the z-component of the vector jV:
The set of the two partial differential equations (32) and (33), together with the boundary conditions (34), is of the so-called ‘‘spatially inhomogeneous reaction–diffusion
Assuming a stress-free surface, rxz = 0, ryz = 0, rzz = 0, relation (26) leads to lV = 0, yielding with lV ¼ Rg T lnðy V =y eq V Þ y V ¼ y eq V
at z ¼ 0:
~y ¼ 1 for f ¼ 0
and
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
(a)
1.010
1.008
3
2
4 5
1.006
yV /y V eq
type” (see e.g. Ref. [16]), and still a topic of ongoing research in applied mathematics; especially the long term asymptotic behaviour is of interest. The problem is, however, more complicated as one would assume at a first glance, since the associated stress field r, reflected by rH, is also not known in advance and evolves together with ~y due to the eigenstrain rate represented by e_ gc ¼ a, see Eq. (5). Since nonlinear terms appear in (32) and (33) as well in the boundary conditions (34), one cannot expect analytical solutions. Therefore, the equations (32) and (33), together with boundary conditions (34), are solved numerically by means of finite difference method on a grid of 101 equidistant points. The time integration is performed by Euler method.
4653
6 7 8 9
1.004
1.002
1.000
1
0.0
0.2
0.4
0.6
0.8
1.0
z/deff
4.2. Simulations of yV and rH profiles
(b)
The simulations are performed for the following values of parameters:
0.0010
1 2
0.0008
3
j
Values, typical for metals, are addressed to 9 v ¼ 3:15Þ. y eq V ¼ 10 ; j ¼ 10 and v ¼ 0:3ð~ The parameter k, which can be denominated as a vacancy intensity factor, relating the diffusion coefficient Deq A and the bulk viscosity factor K, is considered as a free parameter. In the simulations the values k = 1000, 10, 1 and 0.1 are chosen.
~H are chosen as The initial values of ~y and r ~y 0 ¼ 1 and r ~H 0 ¼ 103 , respectively. Four sets of diagrams, Figs. 2–5, are presented according to four distinct values of the vacancy intensity factor k, which obtains k ? 1 for ideal sources and sinks and k ? 0 for no sources and sinks in the bulk. Figs. 2–5 demonstrate the relative vacancy site fraction ~y ¼ y V =y eq V profile across the layer, starting from the surface towards z = deff. A common feature of all the four diagrams is that ~y changes quickly from its initial state (curve label 1), ~y ¼ 1, to its highest values (curve label 2) during the first time increment. Such a supersaturation with vacancies occurs due to the generation of vacancies at sources in the bulk induced by the stress state. This process is studied in detail with much smaller time increments in Fig. 6 for k = 10. One observes that after approximately 20 small time increments, which means approximately for s ¼ 0:4 (a time period much shorter than one distinct time increment in Figs. 2–5) an envelope of the profiles is formed, pointing to a maximum value of ~y ¼ 1:008. Obviously the bulk provides enough sources of vacancies, and the initial equilibrium distribution ~y ¼ 1 is quickly forgotten. However, with an ongoing diffusion process the vacancies move towards the surface, provoke further generation of vacancies causing stress relaxation and thus reduce the supersaturation in the bulk (see Figs. 2–5). After a certain time interval, say 20 time increments, the equilibrium profile ~y ¼ 1 is approached again. This coupled vacancy generation/annihilation and diffusion process in the bulk and the surface acting as ideal sink/ source for vacancies allows for a rather gradual reduction
4 0.0006
5
σH /E
j
6 7
0.0004
8 9
0.0002
0.0000 0.0
0.2
0.4
0.6
0.8
1.0
z/d eff 9 ~H profiles for y eq Fig. 2. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3ð~v ¼ 3:15Þ and k ¼ 1000. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 2.5 107 with respect to the preceding one.
of the eigenstress state, represented by the profile of relative ~H for all four values of k (see Figs. 2–5). hydrostatic stress r The eigenstress state decreases most quickly near the surface. But with ongoing time (curve labels >10) the eigenstress state is significantly relaxed and finally disappears nearly totally across the whole layer thickness. 4.3. The role of nonlinearities in the differential equations One can see immediately from Figs. 2–5 that ~y is close to 1, which allows a replacement of ~y by 1 þ D~y (the symbol D should not be confused with the Laplace Operator D) and a linearization of lnð~y Þ D~y . Then one can rewrite the set of differential equations (32) and (33), ordered with respect to linear and nonlinear terms, as dD~y @ 2 D~y ~H Þ; ¼ kðD~y j~ rH Þ þ nD~y ðD~y ; r ds @f2 d~ rH y eq rH Þ: ¼ V kðD~y j~ ~v ds
ð35Þ ð36Þ
4654
(a)
1.008
2
3
yV /y V eq
1.006
1.003
3
2
4
4 5
5 6
1.002
6 7
1.004
8 9
yV /y V eq
(a)
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
7 8 9
1.001
1.002
1.000
1.000
1 0.0
0.2
0.4
0.6
0.8
1 0.0
1.0
0.2
0.4
0.6
(b)
1
0.0010
1
2
2
3
0.0008
3
0.0008
4 5
5
0.0006
6 7 8
0.0004
6
0.0006
7
σH /E
σH /E
4
9
8
9
0.0004
0.0002
0.0002
0.0000
0.0000
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
9 ~H profiles for y eq Fig. 3. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3ð~v ¼ 3:15Þ and k ¼ 10. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 2.5 107 with respect to the preceding one.
The nonlinear contribution nD~y to the first differential equation (35) follows after rearrangement as ( ) 2 y j~v @D~y @ð~ r2H Þ @ 2 ð~ r2H Þ eq @D~ ~H Þ ¼ y V : þ þ nD~y ðD~y ; r @f 2 @f @f @f2 ð37Þ The corresponding boundary conditions (34) can be reformulated as D~y ¼ 0 for f ¼ 0 and
@D~y ¼ nr~H ¼ 0 for f ¼ 1 @f
0.6
0.8
1.0
z/deff
z/d eff
ð38Þ
with nr~H ¼ j~v
1.0
z/deff
z/deff
(b) 0.0010
0.8
9 ~H profiles for y eq Fig. 4. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3ð~v ¼ 3:15Þ and k ¼ 1. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 5.0 107with respect to the preceding one.
the boundary conditions (39) with simulation results without the nonlinear terms. The results of simulations differ within 0.2% to the simulations including all terms and cannot be distinguished in Figs. 2–6. If the nonlinear terms can be neglected, the differential equation (35) can be reformulated as d @ 2 D~y V~ ~H ¼ D~y þ eq r ; ds yV @f2 d~ rH k eq ¼ y V ðD~y j~ rH Þ ~v ds
ð40Þ ð41Þ
with the boundary conditions @ð~ r2H Þ : @f
ð39Þ
The role of the nonlinear contributions is checked numerically by comparing simulation results including the terms nD~y in the differential equation (35) and nr~H in
D~y ¼ 0 for f ¼ 0
and
@D~y ¼ 0 for f ¼ 1: @f
ð42Þ
Two coupled partial differential equations for the two vari~H ðf; sÞ are to be solved. The following ables D~y ðf; sÞ and r consequences can be drawn:
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
(a)
1.0004 j
2
1.0003
3
4
yV /y V eq
5 6 7
1.0002
8 9 j
1.0001
1.0000
1 0.0
0.2
0.4
0.6
0.8
1.0
z/d eff
(b)
0.0010
1
0.0008
~v 1 @ 2 D~y ¼ eq j~ rH ; yV k @f2 d~ rH y eq y eq rH ¼ k V D~y : þ k V j~ ~v ~v ds
0.0006
ð43Þ ð44Þ
This formulation allows in principle a semianalytic solution of Eq. (43) for D~y via the Green’s function to the operator on the left side and finally the construction of an integro~H . differential equation for r
3 4
σH /E
~H are replaced by D~y and r ~H in the set If D~y and r of differential equations (40) and (41) as well as in the boundary conditions (42), all equations remain unchanged. This means that a solution for the ‘‘tension case” ð~ rH 0 > 0Þ provides automatically the solution for the ‘‘compression case” ð~ rH0 < 0Þ. The quantity ð~v=y eq rH is usually strongly dominant in V Þ~ 9 comparison to D~y , because 1=y eq V is in order of 10 . Since ~ H are smooth, see Figs. 2–5, one the functions D~y and r can neglect D~y on the left side of Eq. (40) yielding by combination of Eqs. (40) and (41)
D~y
2
5 6 7
0.0004
8 9
0.0002
0.0000 0.0
0.2
0.4
0.6
0.8
1.0
z/deff 9 ~H profiles for y eq Fig. 5. Evolution of (a) ~y and (b) r V ¼ 10 ; j ¼ 10; v ¼ 0:3 ð~v ¼ 3:15Þ and k ¼ 0:1. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . . , 21 correspond to values of s increased by Ds = 5.0 108 with respect to the preceding one.
7 6 5
1.006
4 1.004
1.002
0.0
0.2
0.4
0.6
0.8
Let h = h(t) be the temporal change of deff from its initial value and ~h ¼ h=d eff . Two contributions to ~ h can be observed, the first one due to the deposition of atoms at the surface, denoted as ~hv , and the second one only due to the deformation process, described by egc and the corresponding elastic strain state, denoted as ~hd . The rate h_ v follows from the mass balance at the surface, expressed in the reference configuration, see e.g. Ref. [17], Section 4.1, as y ð45Þ h_ v A ¼ jA;Z ð0; tÞ ¼ jV ;Z ð0; tÞ: XA
3
h_ v ¼ XA jV ;z ð0; tÞ
2
If a can be approximated by XA div jv, see Eq. (16), which is acceptable for times s > 0.4, (see Fig. 6), then integration of Eq. (16) and using Eq. (15) provide Z d eff Z d eff @jV ;z dz ¼ XA jV ;z ð0; tÞ ¼ 3 XA ð48Þ e_ gc dz: @z 0 0
1
1.000
4.4. Simulations of the evolution of the layer thickness
Inserting Eq. (17), together with Eqs. (30) and (31), into Eq. (45) provides d~hv y @ eq @~ 2 ~H Þ : ¼ y V ~y ðj~vr ð46Þ ds @f @f ~H in A direct relation between ~hv and h~ rH i, the average of r the interval 0 6 f 6 1, can be derived. Using yA 1 in Eq. (45) one obtains
1.008
yV /y V eq
4655
1.0
z/d eff 9 Fig. 6. Evolution of ~y profile for y eq V ¼ 10 ; j ¼ 10; v ¼ 0:3 ð~v ¼ 3:15Þ and k ¼ 10. The profile 1 corresponds to s = 0, and the further profiles 2, 3, . . ., 21 correspond to values of s increased by Ds = 0.02 with respect to the preceding one.
Then putting Eqs. (47) and (48) together gives Z d eff h_ V ¼ 3 e_ gc dz: 0
ð47Þ
ð49Þ
4656
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657
Differentiation of Eq. (21)3 with respect to time provides 3ð1 vÞr_ H ; 2E
k=0.1
ð50Þ
and its insertion into Eq. (49) yields Z 9ð1 vÞ d eff 9ð1 vÞd eff _ hri: h_ V ¼ r_ H dz ¼ 2E 2E 0
0.0008
0.0006
ð51Þ
Integration of Eq. (51) in time and expression of the variables in dimension-free ones yield
<σH /E>
e_ gc
0.0010
~ ~H 0 Þ: hV ¼ ~vðh~ rH i r ð52Þ ~ The contribution hd due to the deformation process can directly be related to the displacement at the surface, see Eq. (23), as wð0; tÞ ~ : hd ¼ d eff
ð53Þ
Note the negative sign on the right side of Eq. (53), since w is counted as positive in the z-direction. Insertion of egc from Eq. (21) in Eq. (23) yields 3ð1 þ vÞ ~ hd ðtÞ ¼ ð~ rH 0 h~ rH iÞ: ð54Þ 2 It is interesting to note that both contributions to the hd are proportional change of the layer thickness, ~ hV and ~ 0 rH iÞ. Therefore, the reader is provided with a to ð~ rH h~ diagram for h~ rH i for various values of k in Fig. 7. The contributions of both deformation mechanisms are in a fixed hV ¼ ð1 þ vÞ=ð3ð1 vÞÞ. relation, depending only v as ~ hd = ~ 4.5. Limit cases for ideal and no sources and sinks for vacancies in the bulk For the layer with ideal sources and sinks for vacancies in the bulk the vacancy intensity parameter k tends to 1, rH must be limited for any s > 0. Then and the quantity d~ds it follows from Eq. (33) that ln ~y ¼ j~ rH , which is equivalent to lV ¼ 0, see Eq. (10). Its differentiation with respect to s rH 1 d~y provides j~ ¼ d~ds , and Eq. (33) can be rewritten as y ds ~v d~y ¼ kðln ~y j~ rH Þ: j~y y eq V ds
ð55Þ
Insertion of Eq. (55) with ln ~y ¼ j~ rH into Eq. (32) provides ~v ~vln2 ~y d~y @2 ~ ¼ y 1 1þ 2j j~y y eq ds @f2 V 2 ~v @ y y eq @~ 2 @~ þ : ð56Þ ln ~y yV 2j @f @f @f ~v @ @~y ~y ðln2 ~y Þ ¼ 0 for f ¼ d eff : 2j @f @f
Thus, the problem reduces to the only one partial differen~H is tial equation (56) with boundary conditions (57), and r
0.0002
k=10
0.0000
k=1000 2x108
4x108
τ
6x10 8
8x108
9
1x10
Fig. 7. Evolution of the relative hydrostatic stress h~ rH i as function of s for 9 v ¼ 3:15Þ and k ¼ 1000; 10; 1 the parameters y eq V ¼ 10 ; j ¼ 10; v ¼ 0:3 ð~ and 0:1. Note that the curve for k = 1000 de facto coincides for h~ rH i due to ideal sources and sinks.
rH i de facto coincides given by ln ~y ¼ j~ rH , whose average h~ with the case k = 1000 in Fig. 7. For no sources and sinks for vacancies in the bulk of the layer, the vacancy intensity parameter is 0. Since vacancies are neither generated nor annihilated, no change in the volume is possible (see the eigenstrain rate e_ gc ¼ 0 for a ¼ 0). ~H Eq. (33) yields then drH/ds = 0, which means that r remains constant in time. Only the vacancy site fraction profile evolves according to 2 ~2H d~y @2 j~vr j~v @ y y eq @~ 2 @~ ~H þ : ¼ 2 ~y 1 r yV 2 ds @f 2 @f @f @f ð58Þ The boundary conditions remain the same as given by Eq. (34). Again, the problem reduces to the solution of only one partial differential equation (58). 5. Conclusions The results achieved in this paper can be summarized as follows: j
j
j
ð57Þ
0.0004
0
The boundary conditions (34) can be then modified as ~y ¼ 1 for f ¼ 0 and
k=1
A model for simulation of stress relaxation in a layer on a substrate by means of diffusion and generation/annihilation of vacancies has been developed. The model assumes ideal sources and sinks for vacancies on the layer surface, ideal or no sources and sinks for vacancies at the interface and non-ideal sources and sinks for vacancies in the bulk of the layer. The driving force for diffusion and generation/annihilation of vacancies stems from the generalized chemical potentials taking into account the actual stress state according to Ref. [5]. The model leads to a boundary value problem, given by a set of two nonlinear partial differential equations for the vacancy site fraction and the hydrostatic stress. The equations are expressed in dimension-free quantities and solved numerically by finite differences.
J. Svoboda, F.D. Fischer / Acta Materialia 57 (2009) 4649–4657 j
The results of simulations are presented for various relations between the diffusivity and activity of non-ideal sources and sinks for vacancies in the bulk. The analysis is completed by calculations of the change of the layer thickness due to the deformation state and due to the deposition/removing of a film of atoms at the layer surface.
6. Outlook The stress relaxation in surface layers can be performed by several mechanisms [1], and their coupling with the vacancy mechanism seems to be a challenge for the nearest future. One can imagine a coupling of the vacancy mechanism with diffusional and dislocation creep, which could complete the picture on stress relaxation in layers at elevated temperatures. Acknowledgements Financial support by Materials Center Leoben Forschung GmbH (project SP19) and by Research Plan of Institute of Physics of Materials (project CEZ:AV0Z20410507) is gratefully acknowledged. The cooperation of authors is also supported in the frame of COST action P19.
4657
References [1] Freund LB, Suresh S. Thin film materials, stress, defect formation and surface evolution. Cambridge: Cambridge University Press; 2003. [2] Gibbs GB. Philos Mag 1966;13:589. [3] Bhandakkar T, Wei Y, Gao H. Math Mech Solids 2009;14:179. [4] Lindgren LE, Domkin K, Hansson S. Mech Mater 2008;40:907. [5] Svoboda J, Fischer FD, Fratzl P. Acta Mater 2006;54:3043. [6] Svoboda J, Fischer FD, Gamsja¨ger E. Acta Mater 2008;56:351. [7] Suo Z, Kubair DV, Evans AG, Clarke DR, Tolpygo VK. Acta Mater 2003;51:959. [8] Garikipati K, Bassman L, Deal M. J Mech Phys Solids 2001;49:1209. [9] Tsurkov I, Grich WM, Gross TS. Analysis of diffusional stress relaxation in submicron Cu interconnect structures using the model with enhanced vacancy diffusivity in grain boundary region. In: Brebbia CA, editor. High performance structures and materials III. Wessex: WIT Press; 2006. p. 85. [10] Wu CH. J Mech Phys Solids 2001;49:1771. [11] Wu CH. Int J Fract 2007;147:227. [12] Varias AG, Massih AR. J Mech Phys Solids 2002;50:1469. [13] Larche` FC, Voorhees PW. Defect Diffus Forum 1996;129–130:31. [14] Spicer JB, Dikmelik Y. J Appl Phys 2008;104:023528-1-6. [15] Danielewski M, Wierzba B, Bachorczyk-Nagy R, Pietrzyk M. J Phase Equilibr Diff 2006;27:691. [16] Di Francesco M, Fellner K, Markowich PA. Proc R Soc A 2008;464:3273. [17] Fischer FD, Waitz T, Vollath D, Simha NK. Prog Mater Sci 2008;53:481.