{'fraction_non_alphanumeric': 0.034777008533247465, 'fraction_numerical': 0.008211238125905651, 'mean_word_length': 4.618423036333484, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Advances in the numerical treatment of grain-boundary migration: Coupling with mass transport and mechanics

This work is based upon a coupled, lattice-based continuum formulation that was previously applied to problems involving strong coupling between mechanics and mass transport; e.g. diffusional creep and electromigration. Here we discuss an enhancement of this formulation to account for migrating grain boundaries. The level set method is used to model grain-boundary migration in an Eulerian framework where a grain boundary is represented as the zero level set of an evolving higher-dimensional function. This approach can easily be generalized to model other problems involving migrating interfaces; e.g. void evolution and free-surface morphology evolution. The level-set equation is recast in a remarkably simple form which obviates the need for spatial stabilization techniques. This simplified level-set formulation makes use of velocity extension and field re-initialization techniques. In addition, a least-squares smoothing technique is used to compute the local curvature of a grain boundary directly from the level-set field without resorting to higher-order interpolation. A notable feature is that the coupling between mass transport, mechanics and grain-boundary migration is fully accounted for. The complexities associated with this coupling are highlighted and the operator-split algorithm used to solve the coupled equations is described.


Introduction
In a recent paper, Garikipati et al. [1] presented a coupled continuum field formulation for the interaction of electric effects, mechanical response and selfdiffusion. Their approach drew upon earlier work by Larché and Cahn [3,4], Nix and co-workers [5,6], Génin [7], Bower and Freund [8], Xia et al. [9]-to name but a few. In Reference [1], a review of these works and others was presented, computational techniques were developed based on the finite element method, and several initial and boundary value problems involving diffusional creep (Nabarro-Herring and Coble creep) were solved.
The formulation introduced in [1] was extended to interdiffusion, with dopants in silicon as motivation, by Garikipati and Bassman [2]. In the current paper, we present a further extension of that same formulation to account for the interaction of grain-boundary migration with stress-driven self-diffusion and electromigration in polycrystalline solids.
Different strategies used to model grain-boundary motion in a computational setting are described in the literature. Sun and Suo [10,11] developed a twodimensional finite element formulation, based on the idea that the energy dissipated during the motion of the boundary must equal the reduction in the free energy of the system. This formulation was used to study several problems including grain growth in a thin film and the competition between surface grooving and grain-boundary migration. A similar variational formulation was developed by Cocks and Gill [12,13] and used to model the evolution of a large network of grains in two dimensions. Another model was developed by Zhao et al. [14], based on the variational formulation of Reitich and Soner [15] and using the level set method of Osher and Sethian [16].
The level set method is also used in the present paper as an interface-capturing technique. However, the goal of the current work is not merely to simulate grainboundary motion, but to capture fully the interaction between this motion and other microscale phenomena that take place in pure polycrystalline materials, namely stress-mediated self-diffusion and electromigration. This distinguishes the current work from the existing literature.
In the current work, we recast the level-set equation in a simpler form by assuming that the level-set function remains a signed distance to the migrating grain boundary (as with the original level-set equation, the use of an extensional velocity field helps maintain this signed-distance function). Mourad et al. [17] tested the resulting level-set formulation extensively and concluded that it is both accurate and robust despite its remarkable simplicity (for some interface-evolution problems, this approach reduces the original level-set equation, a nonlinear hyperbolic PDE, to an ODE that is almost trivial to solve). They conducted several numerical experiments to assess the ability of the simplified level-set scheme to capture the correct solution, particularly in the presence of discontinuities in the extensional velocity and/or in the gradient of the level-set function. They also examined the convergence properties of the method and its performance in a variety of problems, including curvature flow and problems where the simplified level-set equation takes the form of a Hamilton-Jacobi equation with convex or non-convex Hamiltonian. Discretizations based on structured and unstructured finite-element meshes of bilinear quadrilateral and linear triangular elements were shown to perform equally well. They also found that sufficient accuracy is available through a standard Galerkin formulation without resorting to any stabilization or discontinuity-capturing [18,19] techniques.
In addition, a variant of the simplified level-set formulation mentioned above was employed by Ji et al. [20] for representing the evolution of phase boundaries over unstructured finite-element meshes. Here, we treat a complex coupled problem of which grain-boundary migration is only one facet, in addition to mechanics, self-diffusion and electromigration. We apply this simplified level-set scheme to the problem of grain-boundary migration within this context, and we demonstrate its implementation as an integral part of the wider computational framework used to solve the coupled problem under consideration.
The remainder of this paper is organized as follows. In Section 2, we summarize the formulation for coupled mass transport and mechanics. Then we examine the thermodynamics and kinetics of grain-boundary migration, and show how this phenomenon interacts with mass transport and mechanics in polycrystals. In Section 3, we formulate the grain-boundary migration problem using the level set method and we describe the computational methods used in the implementation of this formulation. Numerical examples are presented in Section 4. A summary is provided and conclusions are drawn in Section 5.

The coupled formulation 2.1 Thermodynamic basis
The thermodynamics is posed in a continuum setting, with motivation provided by atomic processes. A schematic of the lattice is shown in Fig. 1. It shows atoms, vacancies, and a free surface; the latter is a source and sink for vacancies. A grain boundary could serve as a source or sink also. Free surfaces and grain boundaries are treated as regions of finite width, δ s and 2δ gb respectively.

Internal energy density
We consider crystalline materials in which the dangling bonds around a vacancy cause an inward relaxation of the surrounding lattice (see Fig. 1). The resulting vacancy relaxation strain can be expressed as PSfrag replacements t = σn n Figure 1: Schematic rendering of a lattice with atoms and vacancies.
where Ω is the atomic volume, f Ω (with 0 < f < 1) is the volume of a vacancy, C v is the vacancy concentration, C eq v0 is the vacancy concentration at thermodynamic equilibrium under vanishing external stress and 1 is the second-order isotropic tensor. The creep strain resulting from the accumulation or depletion of atoms at a free surface or grain boundary with unit normal, n, can be expressed as A non-phenomenological evolution equation for θ c was derived and discussed in detail by Garikipati et al. [1]. The thermal strain is given by where T is the temperature, T 0 is a reference temperature and α is the linear coefficient of thermal expansion. The elastic strain is obtained by subtracting these inelastic strain contributions from the total strain, ε; i.e.
The stress is obtained from the elastic strain and the (generally anisotropic) fourth-order elasticity tensor, C, as Given a state with local stress, σ, the incremental elastic strain energy density is then given as σ : δε e = (C : ε e ) : δε e , where δε e is the increment in elastic strain. With this background, the incremental internal energy density δe, corresponding to a state {ε, T, C v }, and increments δε and δC v is where η is the entropy density and e f v is the vacancy formation energy. The specific form ofê η (η), the entropic dependence of e, is not important to the development that follows.

External work density
The density of work done by external agents can be expressed as Here, the first term of the right-hand side is the classical stress-power term. The second term accounts for the apparent work performed by the electrostatic potential, ψ, during electromigration; in this phenomenological treatment, q is the apparent charge ascribed to each vacancy. The last term accounts for the work done against the stress when vacancies are created at a free surface or grain boundary. The numerical factor appearing in this term arises from geometrical considerations. Additionally, the requirement that this term be active only at sources and sinks is enforced using the indicator function, χ, defined as

Entropy density
The total entropy density is given by Since the formulation is isothermal, the specific form of the vibrational term, η vib (T ), is unimportant. The second term is the entropy density due to mixing (see Kittel and Kroemer [21] for details), k is the Boltzmann constant, C a is the concentration of atoms, and C s is the lattice site concentration. Assuming that C s remains fixed in any material volume, i.e. δC s = δC a + δC v = 0, the incremental entropy density corresponding to an increment in vacancy concentration, δC v , at given stress and temperature can be expressed as

The Gibbs free energy density
At the given state of stress and temperature, the incremental Gibbs free energy density corresponding to increments δε and δC v is defined as With this, the constitutive relations can be obtained in a systematic fashion as outlined in the following section.
Remark 1 Since free surfaces and grain boundaries are considered, there are accompanying surface and grain boundary energies, γ s and γ gb . These energies are taken to be independent of the strain and vacancy concentration for this formulation, and therefore do not appear in the incremental Gibbs free energy density.

Constitutive relations
Anisotropic elasticity is assumed and the relation between the stress, σ, and the elastic strain, ε e , is given by (5). The relation between the current density, i, and the electric potential, ψ, is given by Ohm's law: where ρ is the electric resistivity.

The chemical potential of vacancies
The vacancy chemical potential is defined in the usual fashion [22]: where µ eq v is a constant reference potential. Applying (13) to (6)(7)(8)(9)(10)(11) gives The coupling with mechanics is evident through the strain-and stress-dependent terms.

The vacancy flux
The vacancy flux is obtained from (14) via the relation where D v is the vacancy diffusivity. This relation has been derived from an atomic basis by Bardeen [23]. The term D v C v /kT is the mobility of a vacancy and −∇µ v , the force acting on it, is the driving force for mass transport.

Governing equations
The constitutive relations established above are incorporated in balance laws for mechanics, electric flow and mass transport, leading to a coupled system of governing differential equations.

Mechanics
Neglecting dynamic effects and body forces, the mechanics problem is governed by the quasistatic equilibrium equation and appropriate boundary conditions: where B is the domain of interest and the boundary subsets ∂B u and ∂B σ have essential and natural boundary conditions specified, respectively. These subsets satisfy ∂B u ∩ ∂B σ = ∅ and ∂B u ∪ ∂B σ = ∂B.

Electric flow
The electric flow problem is governed by the charge conservation equation with Dirichlet and Neumann boundary conditions: where ∂B ψ ∩ ∂B i = ∅ and ∂B ψ ∪ ∂B i = ∂B.

Mass Transport
The mass-transport problem is governed by the continuity equation for vacancies, and appropriate initial and boundary conditions: where ∂B Cv and ∂B j v are the boundary subsets on which concentration and flux boundary conditions are specified, respectively. Here, the boundary subsets satisfy ∂B Cv ∩∂B j v = ∅ and ∂B Cv ∪ ∂B j v = ∂B. The effectiveness of vacancy sources and sinks is characterized by the relaxation time, τ . The equilibrium vacancy concentration, C eq v , is defined by µ v | C eq v = µ eq v in (14).

Grain A Grain B
Grain-boundary region Surface-boundary region 2δ gb δs Figure 2: Variation of e f v and e d v across boundary regions.

Grain-boundary migration
The current location of a migrating grain boundary determines the value of the indicator function, χ(x, t) (see Eq. (8)); i.e. it determines whether vacancy sources/sinks are active and whether creep strain can accumulate at a given point, x ∈ B. Furthermore, information about the location of the boundary is needed to calculate the value of the vacancy formation energy, e f v , and the activation energy for diffusion, e d v , everywhere in the domain of interest. These properties are assumed to vary linearly over the width of a boundary region as shown in Fig. 2.
The foregoing illustrates the influence of grain-boundary migration on mass transport. Due to the tight coupling between mass transport and mechanics, the migration of the grain boundary also affects the stress. In turn, mass transfer across the grain boundary causes one grain to grow at the expense of its neighbor and thus leads to grain-boundary migration.

Thermodynamic driving forces
Generally, interface migration in polycrystals is driven by the accompanying decrease in the free energy of the system. The thermodynamic driving force for such a process, acting on a unit area of the interface, is thus defined as where δG is the increase in the total Gibbs free energy of the system brought about by a motion of the interface, during which the interface sweeps through the volume δV . Neglecting triple junctions, the Gibbs free energy of a polycrystal can be expressed as the sum of two contributions: where g is the Gibbs free energy density as defined in (11) and γ i is the energy per unit area of an interface, Γ i . The term interface is used here to refer to free surfaces as well as grain boundaries, and the summation in (20) is over all such interfaces in the polycrystal.
In situations where the Gibbs free energy density, g, suffers a decrease across a grain boundary, the total free energy of the system can be reduced if the grain with the smaller value of g (evaluated at the grain boundary) were to grow at the expense of its neighbor. Thus, a driving force acts on the grain boundary. For example, during recrystallization, annealed grains grow at the expense of coldworked grains in which large dislocation densities lead to high values of g. The misorientation between two adjacent grains of an elastically anisotropic material subjected to a directional load causes one grain to store a smaller amount of strain energy per unit volume than its neighbor. This leads to strain-induced grain-boundary migration. Electromigration also causes atoms to jump across grain boundaries and thus leads to the migration of these boundaries [24].
Under the current formulation, the driving force discussed above can be characterized as follows: When a single atom hops across the grain boundary, it exchanges positions with a vacancy. Therefore, the volume change associated with this hopping event is δV = (1 − f )Ω. The decrease in the Gibbs free energy of the system is equal to the difference in the chemical potential of atoms across the boundary; i.e. −δG = ∆µ a , or equivalently −δG = −∆µ v . Hence, from Eq. (19), and assuming that the chemical potential gradient across the boundary is essentially linear [25], the driving force, p g , can be expressed as where n is the unit normal to the grain boundary and 2δ gb is its width.
From Eq. (20), it is clear that the total Gibbs free energy of a polycrystal can also be lowered by decreasing the total surface area of the grain boundaries in the system. When a curved interface moves away from its center of curvature, PSfrag replacements γ gb γ gb n R γ gb R Figure 3: Driving force on a cylindrical grain boundary with a radius of curvature, R, due to surface tension, γ gb .
sweeping through an increment of volume δV , the free energy of the system increases by where γ i is the (constant) specific surface energy of the interface and R 1 and R 2 are its principal radii of curvature. This expression can be traced back to Herring [26]. It follows that, in the current two-dimensional formulation, the driving force acting on a unit area of a (cylindrical) grain boundary, due to this effect, can be expressed as where R is the radius of curvature and the negative sign indicates that p γ drives the boundary to migrate toward its center of curvature (see Fig. 3). In a polycrystal comprising only annealed grains, the action of p γ leads to normal grain growth; i.e. the growth of large grains at the expense of smaller ones. By contrast, in the early stages of recrystallization, small annealed grains grow at the expense of the surrounding matrix of cold-worked material. In this case, the boundaries surrounding the annealed grains move away from their respective centers of curvature under the effect of p g .

Kinetic law
Assuming that grain-boundary migration takes place as a result of individual atoms hopping across the boundary, the migration velocity can be expressed in the following form: where v n is in the direction of the local unit normal, n, which is assumed to point away from the boundary's center of curvature and M is the grain-boundary mobility. This classical result, derived by Turnbull [25] using absolute reaction rate theory, holds provided that pΩ ≪ kT , a condition which is always met in grain growth and recrystallization [27]. The results of molecular dynamics simulations of curvature-driven [28] and strain-induced [29] grain-boundary migration in bicrystals agree with Eq. (24), which can be modified as follows, to account for both types of driving forces discussed in Section 2.4.1:

Grain-boundary mobility
In Section 2.4.1, it was established that diffusion of atoms across a grain boundary due to the local gradient in the atomic chemical potential causes the boundary to migrate. The driving force for boundary migration in this case, denoted p g , is given by (21). An expression for the corresponding mobility, M g , can be obtained by stipulating that, in this case, the migration velocity in the direction of the local unit normal, n, should be given by (see Porter and Easterling [30]) where j a and j v are, respectively, the local atomic and vacancy fluxes. Also recall that Ω is the atomic volume and f Ω is the volume of a vacancy. Since atoms and vacancies move by exchanging positions; i.e. j v = −j a , we have Finally, by combining (27), (21) and (15), we obtain Since grain-boundary migration under the sole influence of p g involves transport across the boundary, and since, in the present treatment, mass transport is assumed to take place through the exchange of positions between atoms and vacancies, no lattice sites are transferred across the boundary in this process. By contrast, migration under the influence of p γ consists of the transfer of lattice sites across the grain boundary [27], via the process of atoms detaching from one grain and attaching to the one on the opposite side of the boundary. It is hence unreasonable to expect the mobilities, M g and M γ , of these two distinctly different processes to be the same.
As mentioned earlier, normal grain growth can be attributed to the action of the driving force p γ . The fact that grain growth is a thermally activated process suggests an Arrhenius-type relationship between the mobility, M γ , and the temperature: where the pre-exponential factor, M 0 , and the activation energy for grainboundary migration, e m , are dependent on the misorientation angle and the axis of rotation [31,32].

Computational methods
In this section, we focus on the level-set formulation of the grain-boundary migration problem, and we present a detailed description of the computational techniques used in its implementation. To solve the coupled problem, the level-set formulation is integrated into the computational framework that was first introduced in Reference [1]. This computational framework is based on an operator-split solution scheme and relies on the finite element method to solve the mechanics, mass-transport and electric-flow problems individually. While Reference [1] addresses mainly the physics of stress-driven mass transport in polycrystals, it also contains details of the computational framework used therein. Here, these details are omitted; however, the operator-split solution scheme is outlined briefly in Section 3.2 to show how the level-set formulation is incorporated into the computational framework used.

Level-set formulation
Evolving interfaces can be tracked using the level set method, originally introduced by Osher and Sethian [16]. A comprehensive review of the method and the computational algorithms used in its implementation can be found in [33,34].
Consider an evolving grain boundary, Γ, which divides the domain of interest, B, into two disjoint open subsets, B − and B + . This situation is depicted in Fig. 4. The boundary can be parameterized with the aid of the scalar function φ(x, t), defined on B, provided that the following conditions are satisfied for all t ≥ 0: The term level set refers to a set of points with a fixed value of φ, i.e. an isocontour of φ; the zero level set represents the grain boundary. Accordingly, the unit normal to a given level set can be defined locally as where · denotes the Euclidean norm. This expression can be evaluated at any point on the zero level set to obtain the local unit normal to the boundary, n + 0 . This definition implies that n + 0 always points into the B + region. The evolution of the level-set field is governed by where F n (x, t) is the (scalar) local propagation velocity of the level set passing through point x. To track the motion of the grain boundary, we require that where v n is obtained at any point on the grain boundary from (30). Although this requirement does not place any restrictions on the choice of F n away from Γ, the solution procedure is simplified greatly if F n is an extensional velocity field; i.e. if ∇F n · n + = 0.
This first-order partial differential equation can be solved for F n , at all points x ∈ (B \ Γ), using (34) as a boundary condition. A simpler alternative strategy for constructing the extensional field is discussed in Section 3.1.1 below.
The level-set function, φ, can be initialized as the signed distance from Γ as follows: It is noted that this initial condition satisfies (31). It also implies that initially, Importantly, this desirable mathematical property of the level-set field is preserved when the velocity field is extensional. This can be shown (see [14]) by noting that Combining (38) and (33) and assuming that φ and F n are smooth, we obtain Thus, if F n is extensional (i.e. ∇φ·∇F n = 0) and φ is initially a signed-distance function (i.e. ∇φ = 1 and ∇ ∇φ = 0), we have which implies that (37) holds for all t ≥ 0. Hence, Eq. (33) reduces to which governs the evolution of φ(x, t) from the initial condition (36).
Recall that the mechanics, mass-transport and electric-flow problems are solved using the finite element method, and that part of the coupling between these three sub-problems, on one hand, and the grain-boundary migration sub-problem, on the other, is through the indicator function, χ(x, t). In this setting, the levelset update is computed at each node via a generalized trapezoidal rule: where x A is the position vector of node A and t n+a = at n+1 + (1 − a)t n , with 0 ≤ a ≤ 1. It is noted that when explicit time integration is used (a = 0), the update operation is trivial since v n , and hence F n , are determined from the known solution at t = t n and do not depend on φ(x, t n+1 ). It is also noted that no stabilization is required due to the use of the simplified level-set equation (41) in lieu of (33). The advantages as well as the performance and stability characteristics of the resulting level-set formulation are studied in detail by Mourad et al. [17].
Additionally, χ is defined more precisely as follows: where H(·) is the Heaviside function.

Velocity projection and field re-initialization
The approach adopted here for constructing the extensional level-set propagation velocity field, F n , is based on the notion that F n (x) = v n (p) for any It is clear that the resulting velocity field satisfies (35) when as depicted in Fig. 5. This is not always the case however; for instance, if x ∈ B + (see Fig. 5), then p is such that It follows that, for all x ∈ B + , F n (x) = const., i.e. ∇F n = 0, which clearly satisfies (35) also. The above arguments apply in B − as well.
Due to the accumulation of numerical error, the level-set field may develop perturbations; i.e. ∇φ may deviate from unity in some regions within B. The field must be re-initialized to neutralize these perturbations and retain accuracy by maintaining ∇φ = 1 (see Eqs. (38)(39)(40)). Box 1 shows the algorithm used to implement the velocity projection scheme outlined above and the reinitialization scheme based on Eq. (36).
It is noted that, to preserve the integrity of the solution, the re-initialization procedure should not change the current location of Γ. In other words, the reinitialized field should have the same zero level set as the original (perturbed) field. The re-initialization scheme described here does not satisfy this condition strictly, i.e. it introduces error into the solution. Re-initialization should therefore be used judiciously. This issue is examined in detail by Mourad et al. [17].
An alternative re-initialization procedure designed to minimize this type of error is described in detail in [35,36]. It must be noted however that, in some cases including curvature-driven migration, quadratic convergence in L 2 is achieved with the level-set update formula (41), and importantly, this optimal convergence rate is preserved by the present re-initialization scheme (see [17,37] for details).
Box 1: Velocity projection and level-set field re-initialization algorithm. Here, the global arrays, Φ, F and D hold the nodal values of φ, F n and x − p , respectively. A division of the zero level set which spans a single element is referred to as a segment and is denoted by L 0 . Also, since bilinear shape functions are used, an element can contain only one such segment.
The algorithm in Box 1 was first introduced by Malladi et al. [38] and was previously used by Garikipati and Rao [39]. It is adopted here for its simplicity, despite being relatively expensive-the complexity of the present algorithm is at best O(n L × n np ), where n L is the number of elements intersected by the zero level set and n np is the total number of nodal points in the mesh (see [17]). A more efficient algorithm, such as the fast marching method [33,40] with O(n np log n np ) complexity, could be employed for velocity projection and levelset field re-initialization on Cartesian grids.

Gradient smoothing
The local migration velocity, v n , on the grain boundary is dependent on the local curvature, κ = 1/R, which is defined as follows: From (32) and (37), it is clear that n + = ∇φ, and the curvature can hence be expressed as where n sd = 2 is the number of spatial dimensions. Since the value of φ is updated at the finite element nodes, it is convenient to use the shape functions to evaluate the spatial derivatives in the above expression. However, since bilinear shape functions are used for simplicity and robustness, Eq. (48) cannot be used to evaluate κ directly. To overcome this difficulty, we introduce a 'smoothed' normal vector fieldñ + , weakly related to ∇φ by where w is an arbitrary weighting function. Equivalently, the nodal values of each component,ñ + i , of the smoothed normal vector can be obtained by minimizing the following discretized functional: where n el is the number of elements in the model, B e denotes an element domain, n en is the number of nodes per element and N A is the shape function associated with node A. This leads to a matrix equation of the form Md = f , where d is the global vector containing the nodal values of the componentñ + i . The global mass matrix, M, and right-hand side vector, f , are obtained from the corresponding element arrays via the usual assembly process. The element arrays in this case are given by Finally, the curvature is evaluated as follows: Remark 2 Using this least-squares technique to smooth the stress field leads to a mixed problem in stress and displacement form (see Zienkiewicz and Taylor [41]). By analogy, using this approach to compute the curvature is formally equivalent to a two-field mixed formulation for φ andñ + .
Remark 3 Similar techniques have been used previously to evaluate the curvature of an evolving interface; e.g. see Chessa and Belytschko [42].

Operator-split algorithm
The coupled problem is solved using an operator-split algorithm. The sequence of operations carried out in one time step is shown in Box 2 to illustrate how the level-set formulation is incorporated into this solution scheme.
(2) CONSTRUCT the extensional velocity field, F n , and RE-INITIALIZE the level-set field (see Box 1).

Numerical examples
Numerical results, obtained by solving the coupled initial and boundary value problem, are presented here with the aim of highlighting some of the advantages of the current approach. A 1 µm wide, 2.5 µm long segment of an aluminum interconnect line is modeled. The values of the material parameters used for Al  Fig. 2. † The diffusivity is given by Dv = Dv 0 exp(−e d v /kT ). ‡ A 0 = γ gb M 0 ; also see Eqs. (29) and (30). are given in Table 1. The segment consists of two pure Al crystals separated by a Σ7 tilt grain boundary (38.2 • misorientation about <111>).
The line is assumed to operate at T = 373 K, with a reference temperature, T 0 = 473 K. Rigid passivation material surrounding the line prevents vacancies from crossing the upper and lower boundaries of the domain (see Fig. 6a); i.e. the condition j v ·n = 0 holds at these boundaries. Periodic boundary conditions are imposed on the vacancy concentration at the left and right boundaries of the segment and an electrostatic potential difference, ∆ψ = 0.0021 V, is applied between these two extremities. Vacancies drift along the electric field, E = −∇ψ, pointing to the right.
In the first example, the migrating grain boundary consists, initially, of two straight (planar) sections, which are joined by a circular (cylindrical) section. The straight sections form 45 • angles with the upper and lower boundaries of the domain as shown in Fig. 6a. Since the vacancy formation energy, e f v , is low inside grain-boundary regions where vacancy sources are also present, vacancies accumulate in such regions resulting in high values of the local vacancy concentration, C v . The location of the grain boundary is revealed locally by the maximum-valued contour of C v . Figure 6 shows the evolution of the vacancy concentration contours in the line.    It is clear that, during the first 300 seconds, the central cylindrical section of the boundary migrates to the right; i.e. toward its center of curvature, while the planar sections are relatively less mobile (Fig. 6b). It must be emphasized that the interaction between mechanics, mass transport, electric effects and grain-boundary motion is accounted for, and the driving force for boundary migration, due to stress-driven diffusion and electromigration, is included in the calculations. However, its effect is overshadowed by the dominant driving force due to the curvature of the boundary.
After 900 seconds, the curvature is approximately the same everywhere on the grain boundary (Fig. 6c). At this stage, a steady state prevails and the boundary continues to travel toward the right without undergoing any further changes in shape. It is noted that the level-set calculations remain stable. The boundary remains smooth and does not develop any spurious cusps or ripples. The 45 • equilibrium angles between the grain boundary and the sidewalls are maintained as the solution progresses.
A contour plot of the magnitude of the vacancy flux at t = 300 sec is shown in Fig. 7. It is clear that a strong vacancy flux exists in the vicinity of the grain boundary. A vector plot of the vacancy-flux field in the neighborhood of the grain boundary is shown in Fig. 8 and although the field is complicated in this neighborhood, it can be seen that vacancies tend to drift in the same direction as the migrating boundary, thus preventing the appearance of a 'trail' of vacancies or other oscillations in the numerical solution in the boundary's wake. It is important to note that the flux field evolves continually as the boundary migrates. It is also notable that the formulation captured this aspect of the coupling between mass transport and grain-boundary migration without additional terms being added to the expression of the vacancy flux to account, specifically, for the effect of moving boundaries.
The second example is concerned with the case where a curved grain boundary evolves into a planar configuration to reduce the free energy of the system. Here, the equilibrium angles at the grain boundary-sidewall intersections are set to 90 • and the initial geometry of the boundary is different and less regular than in the first example. It is clear from Fig. 9, which shows the evolution of the vacancy concentration contours in this case, that the boundary flattens as the solution progresses. It is noted that the equilibrium angles are also preserved in this case. It is also noted that in this case, perturbations in the solution lead to the formation of spurious 'shadow' regions ( B + , B − ; see Fig. 5) and re-initialization becomes necessary in these regions to maintain accuracy and to preserve the contact angles.
Although the numerical stability characteristics of the level-set formulationand those of the operator-split algorithm-are not examined in a formal setting, no spatial or temporal oscillations are observed in the numerical solution of the example problems presented, as long as the time-step size is within the CFL limit; i.e. ∆t < h/F n , where h is the mesh parameter. A detailed study of the performance and stability characteristics of the present level-set formulation is  presented by Mourad et al. [17]. Details regarding the advantages, applications and numerical stability characteristics of operator-split schemes can be found in [43,44,45] and references therein.

Summary and Concluding Remarks
In this paper, we present the following contributions.
• A computational formulation capable of capturing the full coupling between grain-boundary motion and other microscale phenomena that take place in pure polycrystalline materials. The coupled continuum formulation presented here was developed to model stress-driven self-diffusion and electromigration in polycrystals while accounting fully for the interaction between these mass transfer processes and the motion of grain boundaries. The formulation accounts for two distinct thermodynamic driving forces acting on a grain boundary; one due to the boundary's own curvature and another engendered by mass transfer across the boundary via stress-driven self-diffusion and electromigration.
• A simplified level-set formulation which does not require spatial stabilization.
The level set method has previously been used to pose grain-boundary migration as a time-dependent field problem governed by a pure advection equation. Standard numerical schemes resort to spatial stabilization techniques (e.g. upwinding schemes, Galerkin/Least-Squares) to attenuate the spurious oscillations known to appear in the numerical solution of equations of this type. Here, on the other hand, the level-set equation is reduced, using the mathematical properties of signed distance functions and extensional velocity fields, to a simpler form which obviates the need for these stabilization techniques. This leads to a remarkably simple explicit scheme for advancing the solution in time. The algorithm used to construct the extensional velocity field-and to simultaneously reinitialize the level-set field-is presented. We also provide the description of an L 2 -projection technique used to compute the curvature of the grain boundary.
The numerical examples presented indicate that the strong coupling in the problem is captured adequately and that the numerical implementation allows the solution of the coupled initial and boundary value problem to be advanced in time in a stable fashion to obtain physically meaningful results.