Modeling Microstructure Evolution in Engineering Materials

ADVANCbS IN APPLIED MECHANICS, VOLUME 76

Modeling Microstructure Evolution in Engineering Materials ALAN C. F. COCKS and SIMON P. A. GILL

and

JINGZHE PAN

I. Introduction.. . . . . . . . . .

..

. ..

11. Microscopic Constitutive Laws . . . A. Diffusional Processes . . . . . . B. Interface Reactions . . . . . . . . C. Migration of Interfaces . . . . . .

.

. . ..

. .... ...... .... . . .... , . .

..

111. Thermodynamic Variational Principle

.

.

.. .. .. ..

. . ....... . . . .

.. .. .. ..

. . . .

.... .... .... ....

...

. . . . 82

. , . . . . . . 84 . . . . . . . . 86 . . . . . . , . 88 . . . . . . . . 89

.. . ...........

.

...

IV. Numerical Models . . . . . . . . . . . . . . . . . . . . . . . . . , . . . . . . A. Grain Boundary Diffusion Controlled Processes . . . . . . . . . . . . . . B. Coupled Grain Boundary and Surface Diffusion . . . . . . . , . . . . . . C. Interface Reactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . D. Coupled Grain Boundary Diffusion and Self-Diffusion . . . . . . . . . . E. Grain Boundary Migration. . . . . . . . . . . . . . . . . . . . . . . . . . F. Effect of Changes in Elastic Stored Energy on Microstructure Evolution V. Rayleigh-Ritz Analyses . . .

.

....

.

.... . .

.

..

..

.....

.

..

90

. 95 . . 96 . . 101

.

. . . .

. . . .

11 1 114 116 136

. . . . . 150

VI. Structure of Constitutive Laws for the Deformation of Engineering Materials . . . . . . . . . . . . . . . . . . . . . . . , . . . . . . . . I54 A. Stage 1 Compaction of Ceramic Components . . . . . . . . . . . . . . . . . I56 VII. Concluding Remarks . . . References . .

....

...

....

.

...

,

.

. .... ......

.

....

.

.

....

..

.

.

.

.

.

....

.

.... ....

. I59

.

. 160

.

...

82

Alan C. F: Cocks et al.

I. Introduction In recent years there has been a growing interest in modeling the microscopic processes that occur within a material when it is processed or suffers inelastic deformation and damage growth under mechanical loading. During processing, it is often important to control certain features of the microstructure. A simple example is the grain size, which can have a significant influence on a range of mechanical and electrical properties (Honeycombe and Bhadeshia, 1995; Spriggs and Dutta, 1974; Heywang, 1971).It is therefore important to develop appropriate kinetic relationships that can predict the evolution of the grain structure. In components processed by a powder route the size and shape of the evolving porosity must also be considered and the influence of this on the grain structure and properties of the finished component evaluated (Du and Cocks, 1992). In alloy systems careful control of the processing conditions are required to ensure that the desired precipitate structure, composition, and distribution are obtained Lo acliicve optimum mechanical or electrical performance (Matan el al., 1998). Similarly, in thin-film devices the exact composition and geometry of the film or islands needs to be controlled to obtain the desired optical or electrical properties (Howes and Morgan, 1985). The microstructure of a material can continue to evolve in the design environment. Dislocation and precipitate structures can change; voids and cracks can nucleate and grow, leading eventually to failure of the component (Cocks and Ashby, 1982). Similarly, the structure of thin films can evolve further: For example, they can break up into islands or fail by the growth and migration of voids (Amerasekera and Nam, 1997; Suo, 1996, and references therein). Before developing modeling strategies to deal with the above situations, it is important to ask what the purpose of the modeling exercise is and to determine what information is available about the fundamental properties of the system being evaluated. For example: Is the aim to exactly model the way in which a given system evolves under a prescribed loading history and in a given environment; or are only the major features of the evolution process required, together with an identification of the dominant kinetic processes; or is the aim to develop a macroscopic constitutive law for the material response in terms of a limited number of state variables? It is also important to determine what is considered to be the fundamental starting point for the development of the models. The terms microscopic andfundamental properties mean different things to different research communities. In ab initio simulations the fundamental properties are the quantum states; in atomic simulations they are the interatomic potentials; in microplasticity models they are the properties of individual dislocations; while in continuum plasticity

Modeling Microstructure Evolution

83

and lunetic models we think in terms of local continuum properties such as yield strength and diffusivities. The more microscopic the starting point the smaller the scale of problem which can be sensibly considered. We also need to consider the accuracy within which the fundamental properties are known. It is not sensible to undertake a large-scale simulation if an accurate measure of the fundamental properties is not available. Here we start from a description at the continuum micromechanical level and describe a general strategy for modeling microstructure evolution which can be adapted to deal with any of the situations identified above. The approach is based on a fundamental variational principle which considers the competition between all possible dissipative processes that are driven by a number of thermodynamic driving forces. Basic microscopic constitutive laws for the different kinetic processes are presented in Section I1 and the variational principle is derived in Section 111. In Section IV we show how numerical schemes can be developed from the variational principle. A number of examples are given which involve a range of kinetic and dissipative processes, such as grain boundary, lattice, and surface diffusion, and interface migration, driven by changes in the total potential energy of the system and changes in surface, grain boundary, and interface energy. By examining the form of the variational principle, it is possible to identify the dissipative processes that will dominate in a given situation. Also the approach allows for a range of different modeling and numerical strategies. The utility of different approaches is discussed for a range of practical physical situations. Even by adopting a reasonably large intrinsic length scale for our starting point, i.e., starting from a micromechanical continuum description, the computational requirement of problems of practical interest can rapidly grow beyond the capabilities of currently available computers. Two different strategies can be adopted to overcome this problem. In the first of these a number of representative fundamental problems can be identified and detailed simulations conducted. Thermodynamically consistent constitutive evolution laws can then be obtained by identifying major features of the evolution process (e.g., in grain growth studies one might examine features of the distribution of grain sizes within the analyzed region) and feeding these characteristics back into the variational principle. This provides evolution laws in terms of a reduced number of degrees of freedom. Large-scale problems can then be analyzed by treating these new relationships as the fundamental relationships for this new situation. The second approach differs in detail but not in philosophy. From the detailed analysis of a fundamental problem, major geometric features of the evolving microstructure can be identified. It is then possible to identify a very coarse description of the microstructure which retains these essential global features. Again, use of this coarse description

84

Alan C. F: Cocks et al.

for the local problem readily allows much larger scale problems to be analyzed. Both these approaches naturally provide a mechanism of smoothing information up the hierarchy of length scales. They also prompt the investigator to think about the physical situation under examination, to consider the purpose of the analysis being conducted and to identify the overall objectives. These procedures are fully described in Sections IV and V. The variational principle employed in this paper also allows general features of the structure of constitutive laws for the macroscopic response to be identified which consider the competition between the different mechanisms. By necessity, these constitutive laws need to be expressed in terms of a limited number of state variables and these procedures can be thought of as an extension of those described above for analyzing large-scale problems. We examine these general features in Section VI and describe a range of constitutive laws which have been developed for creeping and sintering materials. 11. Microscopic Constitutive Laws

In this section we develop a strategy for analyzing a wide range of physical processes. We seek to describe the micromechanical physical processes in terms of continuum-type properties. Typical examples of the range of kinetic and dissipative processes we are interested in are shown diagrammatically in Figure I and listed in Table I . These principally involve the diffusional rearrangement of material within a body and/or the motion of internal boundaries and interfaces. The thermodynamic driving forces for each of these processes are identified below. There is no unique starting point for the analysis of physical situations in which these mechanisms operate and which fully take into account the interactions between the different mechanisms and driving forces, but any procedure must be thermodynamically consistent and allow the relative importance of all the competing mechanisms to be identified. Our objective is to determine the macroscopic response of a body. We consider an element of material within the body as a thermodynamic open system, which can exchange material with surrounding elements and whose state can change as a result of these exchanges. The macroscopic response can then by obtained by integrating over all the elements. Before doing this we need more information about the constitutive behavior at the microlevel. At this stage we wish to express the material response in as general a way as possible. To clarify our exposition of the method of analysis proposed here, however, it proves useful to limit the range of kinetic processes, number of diffusing species, and different types of physical problems in which we are interested. In particular, we assume that there

Modeling Microstructure Evolution

Ti

t t t t t

85

r-----------

1 1 1 1 1 FIG.I . General situation considered iii this paper. The microstructure can evolve by a number of kinetic/dissipative processes: grain boundary. lattice. or surface diffusion: the migration of interfaces: evaporation and condensation of material: or the rate of operation of any interface reactions.

is always an equilibrium concentration of vacancies in the matrix, so that the diffusivities are not a function of state, and we limit our attention to situations in which diffusion of a single element determines the response. We do, however, allow for phase changes, which can be accompanied by a change in volume of the material. It is possible to extend the approach to a wider range of situations and this is discussed more fully in subsequent sections.

TABLE1 F U L LR A N G EOF KINETIC/DISSIPATIVE PROCESSES AN11 THERMODYNAMIC DRIVING FORCESC O N S I D E R E I ) IN THIS PAPER KineticlDissipative Mechanisms ~~

Origin of Thermodynamic Driving Forces

~

Lattice diffusion Grain boundary dilfusion Surface diffusion Interface reactions Grain boundary migration

Potential energy of applied load Stored strain energy Surface energy Grain boundary energy

86

Alan C. F: Cocks et al.

It proves convenient in this section, and in the derivation of the variational principle in Section 111, to express the material response in terms of the transport of atoms within the body. When considering practical problems in Sections IV to VI, it is generally more convenient to express the behavior in terms of continuum quantities. In particular, in the absence of phase changes, it proves more convenient to consider the transport of volumes of rnaterial from one location to another. The normalizations and groupings of material parameters adopted in this section are dictated by this practical utilization of the constitutive laws.

A. DIFFUSIONAL PROCESSES The constitutive behavior of all the diffusional processes identified in Table 1, namely lattice, grain boundary, and surface diffusion, can be represented in the same general form. For lattice diffusion we consider the flux of atoms across unit area within a body, J l , where i has a value in the range I to 3 and refers to a right-handed coordinate system. According to Fick’s first law (Shewmon, 19631,

where w is the excess chemical potential of an atom at position x;, Dl = D / Q / k T is an effective diffusivity, with D/ representing the lattice diffusivity, which, in general, is a function of the state of the material, Q is the atomic volume, k is Boltzmann’s constant, and T is the absolute temperature. The reference state for definition of the chemical potential is a flat stress-free surface. The flux is driven by the potential gradient s,.We will see later that it is instructive to express this relationship in potential form, such that

where $1

1 Dl

= 5 $SI.Sl.

The rate of energy dissipation per unit volume for material rearrangement due to lattice diffusion 21 can be expressed in terms of the potential &:

(2.3)

87

Modeling Microstructure Evolution is the dual potential to q5! and si

= --.a @I

a J;

When analyzing situations where material rearrangement occurs by diffusion in an interface, such as along a grain boundary, along a free surface, or along the interface between two different solid phases, we consider the material to flow through a thin interface layer of thickness &in.It is then convenient to express the diffusional constitutive laws in terms of the atomic flux across unit length of interface. Then (2.4) with

The effective diffusivity DInis now related to the interface diffusivity D,, through the relationship V,,= D1nSlnR/kT. In each of these expressions we have represented the response in terms of a local coordinate system, x;, in the plane of the interface, such that the subscript a has a value in the range 1 to 2. Throughout this paper we use Greek subscripts to refer to a coordinate system within the plane of an interface. A repeating suffix then implies summation over two components. Arabic symbols refer to a general three-dimensional coordinate system and in this instance a repeating suffix implies summation over all three components. If n is the outward normal to the material on a given side of an interface, the chemical potential of an atom of that phase in the interface is

( - a , J n l n / - ? 4 n ( K l + K 2 ) +UP)" = (- fin(K1 + K 2 ) -k "c)",

P =

(2.5)

air

where a,, is the stress normal to the boundary, vIl is the free energy per unit area of interface, K I and K? are the principal curvatures of the interface, which are positive for concave surfaces with respect to the outward normal, and u p is the strain energy density. For a free surface, a,]is simply - p , where p is the pressure exerted by the vapor phase on the surface, which is constant for a given continuous surface. For a grain boundary, the phases are the same on each side of the boundary, with the interface curvatures equal in magnitude but opposite in sign. The net flux of material in the boundary is then only dependent on the stress normal to the boundary and the strain energy density, i.e., p = (-a,l u,)".

+

Alan C. F: Cocks et al.

88

When considering energy dissipation it is convenient to consider the energy dissipated per unit area of interface d i n . Now d i n = @in

+ +in,

where

is the dual potential to @in and s, =

a +in a J,

-*

B. INTERFACE REACTIONS In the above expressions for diffusion in an interface, we have implicitly assumed that these interfaces are perfect sources and sinks for the diffusing species; i.e., all the work done during the diffusion process is available to drive the diffusive flux of material. In practice, some of this work is used to drive the sources and sinks. Following Ashby (1969), Cocks (1992) took this into account by modifying eq. (2.4). He identified part of the total chemical potential with the diffusion process and part with the deposition process, so that, for diffusion in an interface,

Interface reactions can also be important when material diffuses through the grains. The effect is then to change the chemical potential boundary conditions at the interfaces where material is either removed or deposited. We will refer to F, as an interface reaction potential. If atoms are added to an interface of unit area at a rate N,,, the rate of energy dissipation is dr = -N,,wr. Cocks (1992) proposed that p r depends on the rate of plating:

where 4,. is a scalar positive function of p r . This ensures that the rate of energy dissipation is positive. Here we will assume a simple power law function of the form

a,.

Modeling Microstructure Evolution

89

where fi" and are suitable reference quantities. Interface reactions determine the creep response of fine-grained materials at low stresses. Chen and Xue (1990) obtained creep exponents of approximately 2 when fitting a Norton law to creep data for a number of fine-grained engineering ceramics. This observation implies that the interface reaction index m is 2. A value of m = 2 is also consistent with the mechanistic models of Burton (1972) and Arzt et al. (1983). The rate of energy dissipation per unit area is now given by dr = @in

+ + + +in

@r

@rq

where

and

OF INTERFACES c . MIGRATION

The microstructure of a material can change by material jumping across an interface, resulting in migration of the interface. Consider the situation where an interface exists between two phases, which we designate as A and B . We define the normal n to the interface such that it points away from phase A and into phase B . If Mi,, is the mobility of the interface, then the velocity in the direction of n is

(2.10) where f is the thermodynamic driving force, which is simply the change in free energy that results when an atom is transported across the interface from B to A . Changes in interface energy, elastic stored energy and chemical energy, all contribute to the thermodynamic driving force. Then

where AM, is the change in elastic stored energy experienced by the element of material and AwCis the change in free energy associated with the phase change at constant temperature and pressure, 'i2 is the atomic volume of phase A , and A f i is the difference in atomic volume between phases B and A (the elastic deformation required to accommodate this volume change contributes to the change in elastic stored energy A u e ) .

Alan C. E Cocks et al.

90

For situations where the two phases are the same, but of different crystallographic orientation separated by a grain boundary, and there are no internal stresses,

f = (KI

fK2)YbQr

where y ~ is , the grain boundary energy per unit area. When modeling phase transformations in materials, all three contributions to the driving force must be considered. Equations (2.10) and (2.11) are also valid for the analysis of material rearrangement by evaporation and condensation. Material rearrangement in the vapor phase, which we take as phase B , is then much quicker than diffusional processes in the solid phase, so that a constant chemical composition and free energy can be essentially maintained throughout the vapor. For a deposition process at constant pressure and temperature, Au, is simply the elastic strain energy in the solid at the surface and, provided there are no compositional or phase variations, Auc is constant. We can express the governing equation for interface migration employing the same potential structure as for diffusional rearrangement, i.e., (2.12) and the rate of energy dissipation per unit area is d1n

where the dual potential

$,n

= 4m

+

$,,I

,

is now given by (2.13)

111. Thermodynamic Variational Principle

In the previous section we identified a number of kinetic and dissipative processes for the evolution of the microstructure. For a given mechanism and driving force, it is possible to develop a number of approaches for analyzing the response of a system, which will all produce essentially the same results. Suitable approximate procedures can also be readily identified which capture the major features of the evolution process and which are thermodynamically consistent. It is then a matter of personal choice as to which procedure should be adopted for a given problem. In situations where a range of mechanisms and driving forces are involved, greater care is required in the development of analysis procedures and the

Modeling Microstructure Evolution

91

status and appropriateness of any approximate procedures need to be carefully assessed. Sun et al. (1996) have demonstrated that a common assumption that the microstructure evolves in a direction which maximizes the change in Gibbs free energy is generally incorrect. Both the thermodynamics and the kinetics of the process determine the direction of microstructure evolution. We illustrate this result in Section 1V.A where we examine material behavior when grain boundary diffusion is the dominant mechanism of material transport. It is evident from the microscopic constitutive laws presented in the previous section that a given change in microstructure, e.g., a change in surface profile, can result from a number of kinetic processes. It is obviously important to identify the appropriate kinetic process with this change and to determine under what conditions a given mechanism dominates. In this section we derive a general variational principle which naturally takes into account the competition between the different kinetic processes and which can be used to aid the development of computational simulation procedures and approximate methods of analysis. The full utility of this variational principle is explored in subsequent sections. Consider an isothermal body of volume V which is subjected to constant tractions T, over part of its surface rr and displacement rates over the remainder rl,. We assume that the microstructure has evolved to a particular state and we wish to determine the rate of change of state and the inelastic deformation rate of the system. We isolate an elemental volume of material, d V , which consists of a single phase and which we consider to be a thermodynamic open system. For simplicity, we only consider the situation where there is one diffusing species, but it can adopt a number of phases. We can identify two types of elements: elements that lie wholly within a grain and elements that have part of their boundary coincident with a free surface or interface. The specific internal energy of an elemental volume is e = e(ei,, s ) ,

(3.1)

where e;,j is the total strain, measured with respect to a suitable reference configuration, and s is the specific entropy. The Gibbs free energy of the entire body is then given by

where ui are the surface displacements, p is the density of the material, and T is the temperature of the body. By considering the different contributions to the free energy, we can relate the rate of evolution of the microstructure to the rate

Alan C. I? Cocks et al.

92

of change of the Gibbs free energy. This is the purpose of the remainder of this section. For any arbitrary change of state accompanied by a change in the number of atoms d N of an element, the internal energy changes by an amount

dEdV =depdV+pdN= = (gij dE;j

+ T ds)pdV +p d N .

In the above equation we identify r ~ ; ,and p with equilibrium stress and chemical potential fields, oi?;.and /A*, in the body which satisfy all the imposed boundary conditions and are continuous functions which are differentiable to second order, and d e ; j , d s , and d N with d t , Y d t , and N " d t where d t is an arbitrary increment of time and the subscript c indicates that the fields k:), ?', and N" are compatible but not necessarily the exact fields. Material can enter an element by diffusing through the grain, along an interface layer or by jumping across the interface. We divide N into three components N i ' d V , N/il d A , and &; d A , representing these three respective contributions. Dividing by dt and integrating over the volume and total interfacial area Ail, gives

where E is the rate of change of internal energy associated with the assumed fields and p; is that portion of p* required to drive any interface reactions. The atomic flow rate into an element of material is related to the flux in the grain

a J, -+N, ax;

=o,

and the rate of plating on the interface due to the diffusional transport of material in the interface is related to the flux through the interface layer

a

JLY

-

ax:,

+ $'il

= 0.

(3.3)

93

Modeling Microstructure Evolution

Substituting these continuity conditions into (3.2)and making use of the principle of virtual work and the divergence theorem gives

s,

T,*ui d r

EdV =

+

J:

+T

s,

p i ' dV

% dV + S,,n(p* p : ) N i d A dx, -

(3.4)

- ip*J:n,dT.

where r is the total surface area with outward normal n and J:nj is the flux of material out of the body. In (3.4) Ain is the total interfacial area for all phases. A grain boundary is an interface to two phases and it contributes twice to the above equation. If we now identify the stress and chemical potential fields with the exact fields, q,,p. and p,., and the assumed compatible field with arbitrary increments d u f , d i e , d J f , dN;,, d J i , d N h , which satisfy any prescribed boundary conditions and which result in a small change in the rate of change of specific internal energy 2,rearranging (3.4) after noting the constitutive laws of Section I1 gives

=-

JI.

d$/ dV -

S,,.

d@md A -

S,,,

d$m d A -

S,,,

d$r d A = -d'JJ,

where G is the rate of change of Gibbs free energy of the system and \I, =

$1

dV

+ S,

,n

$rn

dA

+ S,, @indA+

S,). $r

dA.

(3.5)

Thus the exact flux and velocity field provide a stationary value of the functional

n,. = w + G .

(3.6)

It can further be shown that for any class of assumed compatible flux and velocity fields there is only one stationary value and this represents a minimum value

Alan C. E Cocks et al.

94

of n,. In situations where our range of assumed fields is limited, we assume that the most appropriate field, i.e., that which is closest to the exact field, is that which provides the minimum value of n, within the assumed set of fields. Suo (1996) has demonstrated how this variational principle can be employed to analyze a wide range of problems, in which the behavior is expressed in terms of a limited number of geometric parameters. In this paper we mainly concentrate on the development of numerical procedures for large-scale simulations, but in the process also examine how simple geometric descriptions can be employed to provide information about the way in which the microstructure evolves. In order to use this variational principle, we do not need to specify the exact form of the constitutive laws described in Section I1 nor determine all the contributions to the chemical potential. These are provided by the variational principle and the exact solution also satisfies all the necessary internal equilibrium requirements. We discuss these conditions more fully in Section IV when we consider the macroscopic response in more detail. An alternative stresdchemical potential-based variational principle can be obtained by identifying u f , iC, J:, N;, J i , and NFn with the exact field u i , i,J j , N m , J a , and &in and o;, w*, and with arbitrary increments doij, d w , and d w r . We now find that the exact stress and chemical potential field is that which minimizes the functional

n,=a+G, where =

S,

4, d v +

(3.7)

lin + S,'" 4m

d~

4in d~

+

Lin

4r

d ~ .

(3.8)

We can further show that the rate of energy dissipation b is given by

D

= -G = @

+ 9.

(3.9)

In practice, it is much easier to base numerical schemes on the flux-based variational principle represented by (3.5) and (3.6) than the stredchemical potentialbased variational principle of (3.7) and (3.8). In the following section we therefore base our numerical schemes on this first variational principle. In many practical situations it proves more appropriate to express the variational principle in dimensionless form. Rather than provide a general normalization here, it proves more convenient to identify appropriate dimensionless quantities when analyzing a given class of practical problems. We return to the variational principle of (3.7) and (3.8) in Section VI where we examine the structure of constitutive laws for the inelastic deformation of engineering materials.

Modeling Microstructure Evolution

95

In this section we have considered a reasonably wide range of mechanisms which encompass the full range of situations considered in this paper. It would have been possible to consider a much wider range of situations, involving the diffusion of a number of elements, which interact to form a number of different phases, and other dissipative processes, such as viscous flow or creep of the grains. It can readily be demonstrated that the general forms of the variational principles represented by (3.5) to (3.8) are still valid for this more general problem and, provided information about the effect of changes of composition on the Gibbs free energy is available, it is possible to extend the procedures described here to this more general class of problem.

IV. Numerical Models In this section we describe how the kinematic variational principle presented in the previous section can be used to aid the development of numerical procedures for modeling microstructure evolution in engineering materials. Throughout this section we limit our attention to the diffusion of a single component and do not consider phase transformations. It then proves convenient to consider the volumetric flux rate of material across unit area, ji = QJi, or unit length of interface, j , = nJ,, rather than atomic fluxes, and the rate of plating of material onto interfaces, u,, = n&,,,rather than the rate of addition of atoms. We start by considering the response of a network of rigid grains in which material rearrangement occurs by grain boundary diffusion, driven by changes in the potential energy of the applied loads. We then gradually increase the number of kinetic processes and thermodynamic driving forces and examine the interaction between these for different classes of problems. The examples described in the following subsections are not exhaustive. They have been chosen to illustrate the proposed methods, to demonstrate that they produce accurate results by comparing them with either available analytical results or physical observations, and to show how the simulations can be evaluated using the variational principle to provide simple macroscopic models of the way in wluch the microstructure evolves. Also, equal emphasis is not given to all mechanisms and the reader is referred to the cited references for more detailed descriptions of the techniques and a more comprehensive evaluation of the results of the simulations. The variational principle described in Section I11 allows the rate of evolution of the microstructure to be determined. In the simulations described in this section we update the microstructure incrementally using a simple explicit Euler scheme by multiplying the velocities by a suitable time step.

96

Alan C. E Cocks et al.

FIG. 2 . A two-dimensional network of grains taken from a micrograph of a ceramic component. Reprinted from Pan and Cocks (1993a). with permission of Elsevier Science.

A. GRAINBOUNDARY DIFFUSION CONTROLLED PROCESSES For situations in which grain boundary diffusion is the dominant kinetic process and the thermodynamic driving force arises from changes in potential energy of the applied load, the variational functional of (3.6) becomes

where Ab is the total area of the grain boundary, L b is the total length of the boundary that intersects a free surface, and m is the unit normal along the the line of the boundary where it meets a free surface, such that it points away from the body. In this section we consider the class of problem represented by Figure 2 in which each grain is bounded by straight boundary facets and all boundaries to the element are axes of symmetry. The last term of eq. (4.1) is then 0. It proves convenient to express eq. (4.1) in nondimensionalized form. If d is a characteristic physical length scale for the problem under consideration, such as the grain size, and a0 is a suitable reference stress, we can identify a characteristic strain rate

97

Modeling Microstructure Evolution

FIG.3. The degrees of freedom and constraints associated with a grain boundary element.

All lengths I , stresses CT, velocities u , grain boundary fluxes j , and time t , can then be normalized such that -

I=-

1 d’

(T=-

0

07,

,

u u=hod ’

;

J = -

J

iod’ ’

i = tho.

(4.3)

Then

and the exact solution occurs when

an,. = 0.

(4.5)

Cocks (19901, Cocks and Searle (1991), and Pan and Cocks (1993a) have described how numerical schemes can be developed from (4.4) and (4.5)when the grains are hexagonal, and Pan and Cocks (1993b) have generalized these schemes to deal with arbitrary two-dimensional shapes of grains. The flux at any point of the boundary can be expressed in terms of the velocities of the centers of the grains on either side of a grain boundary (two translational velocities, i f and , one rotational velocity, w , for each grain in two dimensions) and the flux at the center of the boundary, as shown in Figure 3, by using the continuity condition of eq. (3.3). The contribution to the variational functional for a given grain boundary facet can therefore be presented in terms of seven discrete quantities, and the first term on the right of eq. (4.4) can be determined by summing the contributions

Alan C. E Cocks et al.

98 from each individual facet:

where [u]is the matrix of grain velocities and fluxes at the center of each facet, [ B b ] is a matrix relating the flux to the degrees of freedom, and [ K b ] is the assembled viscosity matrix for grain boundary diffusion. In the derivation of (4.6) we solved a series of local problems and did not examine how these local problems interact with each other. Matter conservation is only guaranteed along a boundary and not at its ends where it is connected to other boundaries. We must therefore further constrain the range of admissible fields by requiring that the flux of material into a triple-point be 0; i.e., at each triple-point, boundaries

where the summation is over all intersecting boundaries and now m, are the components of the unit normal along the line of the intersecting boundary pointing away from the triple-point. We account for these constraints by introducing a series of Lagrange multipliers, one per triple-point, so that the variational functional of (4.4) becomes

points

where the outer summation is over all triple-points and h are the Lagrange multipliers, which physically represent the values of -ii/fi at the triple-point. Thus the chemical potential is continuous throughout the network of grains. The advantage of enforcing conservation of matter using Lagrange multipliers in this way is that the resulting viscositykonstraint matrix can be organized into a narrowly banded matrix, for which there are many efficient solvers. An alternative approach would be to directly relate the fluxes at the triple-points to each other. This would be much more difficult to formulate and would result in a dense viscosity matrix, although there would be fewer variables for which to solve. We can express the functional of (4.7) in discretized form: fic

1 = j[UI'[Kbl[u]

- [UlT[FI

+ [AlT[Cbl[UI,

(4.8)

where IF] is the matrix of applied forces, [ c b ] is a constraint matrix relating the fluxes at the triple-points to the velocity and flux degrees of freedom [u] through matter conservation, and [A] is the matrix of Lagrange multipliers. Differentiating

99

Modeling Microstructure Evolution

FIG.4. Superplastic deformation of an m a y of grains. Reprinted from Pan and Cocks (1993a). with permission of Elsevier Science.

(4.8) with respect to [u]and [A]to find the stationary value of I?,

gives

This equation represents a set of linear simultaneous equations wluch can be solved using standard matrix methods, providing the exact values of [u]and [A] for the posed problem. Care needs to be taken, however, in the choice of method, because the leading matrix is not positive definite. The geometry can be updated by multiplying the calculated velocities by an appropriate time step and adjusting the positions and relative orientations of the grains. In the process, the lengths of the boundary facets will change. As the microstructure is allowed to evolve, grains can change neighbors and facets, and eventually grains can disappear. Methods of dealing with this are fully described by Pan and Cocks (1993b). These topological changes are similar to those that occur during normal grain growth by grain boundary migration. These are described more fully in Section IV.E.l. Pan and Cocks (1993b) have used this method to model superplastic deformation in engineering ceramics. A simulation starting from a grain structure taken from a micrograph for alumina is shown in Figure 4. All boundaries were assumed to be axes of symmetry and the simulation was continued until the structure experienced 100% strain. Near the boundaries the grains experience significant elongation as a result of the constraint imposed by the boundary conditions, but away from the boundaries, where the grains are able

100

Alan C. E Cocks et a1

to slide and rotate over each other, an equiaxed grain structure is maintained. As the microstructure evolves, 8% of the grains are lost and there are many grainswitching events. A small cluster of grains is highlighted in Figure 4. These have changed neighbors and become strung out in the axial direction as the body has been deformed. The results of this simulation are consistent with the major observations of the way in which the grain structure evolves during superplastic deformation of ceramic materials. The grain-switching events do not, however, have a significant effect on the overall deformation process, contrary to the model proposed by Ashby and Verrall (1973), and the strain rate is well represented by a simple Coble creep analysis for an array of regular hexagonal grains with the same mean grain size (Cocks and Searle, 1990). Simulations of this type provide valuable information about the details of the way in which the microstructure evolves and the geometric and material features that most strongly influence the material response. If the major objective, however, is to determine the overall, or average, response, then analyses based on simple geometric models, which can often be performed analytically, generally provide an adequate description of the material response. In the above description we limited our attention to situations in which the grains are formed by a small number of straight facets. In practice, the boundaries might be curved. The above procedures can be readily extended to this situation by discretizing a curved boundary by a number of straight elements. We follow exactly the same procedure as used above, relating the flux along each element to the velocities of the adjacent grains and the flux at the center of the element and introducing Lagrange multipliers at each nodal point formed where two elements meet to ensure continuity of flux (see Pan et al., 1997, for further details). The set of simultaneous equations represented by (4.9) was obtained by considering both the dominant kinetic process and the thermodynamic driving force. We obtain similar relationships for all the situations considered in this paper. In Section I11 it was noted that a common assumption employed when modeling microstructure evolution is that the microstructure evolves in the direction that maximizes the rate of change of Gibbs free energy. It is instructive to examine how the governing equations obtained from such a model would differ from those obtained from the variational principle. Consider the situation where the evolving geometry can be represented in terms of a number of state variables, x l , such that G = G(x,) and

where Fj are the thermodynamic forces associated with the state variables x i . The rate of change of Gibbs free energy is maximized if Xi is normal to the surface of

101

Modeling Microstructure Evolution constant G in state variable space, i.e., if .ii

(4.10)

=Me,

where M is a constant, which is determined by the lunetic processes. Equation (4.10) implies that the viscosity matrix is simply the identity matrix multiplied by a kinetic constant, allowing the rate of evolution of each state variable to be determined independently of each other. This is obviously an attractive feature of this approach, as it permits simple and efficient numerical schemes to be developed. Problems can, however, be encountered with the implementation of internal compatibility requirements, and often ad hoc adjustments are required to ensure continuity of the microstructure. Numerical schemes based on (4.10) only provide an approximate model of microstructure evolution. Schemes based on the variational principle of (3.6) provide the most appropriate model for a given discretization of the real problem. So far we have only considered situations in which there is a single kinetic process. In the following sections we examine situations in which more than one dissipative process contributes to the way in which the microstructure evolves. It is not immediately evident what the single most appropriate kinetic constant is in eq. (4.10) when assuming that the microstructure evolves in the direction of maximum G.The variational principle of (3.6),however, naturally takes into account the full interaction between the different lunetic processes. We therefore employ the variational approach when modeling these more complex situations.

B.

COUPLED

GRAINBOUNDARY AND

SURFACE

DIFFUSION

In this section we examine situations in which both surface and grain bound-

ary diffusion contribute to the evolution of the microstructure. We also include

the influence of changes of surface and grain boundary energy on the evolution process. Employing the normalizations of (4.3), the variational functional of (3.6) becomes

(4.11)

102

Alan C. F Cocks et al.

FIG. 5. The classical model of Hull-Rimmer void growth. A represenlative cylindrical unit ccll of radius I centered on a grain boundary pore of radius r/, is isolated from the body and assumed lo experience the same stress cr that is applied to the body.

where

- the total normalized grain boundary and surface areas, The areas A h and A,y are which evolve at rates A and A , F . The associated changes in total grain boundary and surface energy provide additional contributions to the change in Gibbs free energy of the body.

1. Fast Su$ace Diffusion

It proves instructive to first consider some simple geometric situations and extreme material parameters using (4.11) before examining the full interaction between surface and grain boundary diffusion. If fi,. is large, i.e., if the surface diffusivity is much faster than the grain boundary diffusivity, the second term of eq. (4.11) only makes a significant contribution to the functional if the surface flux is much faster than the flux along the grain boundaries. If they are comparable, then the influence of surface diffusion on the global evolution process is small. It proves useful to discuss this point further by examining a simple problem. Figure 5 represents the classical model of Hull-Rimmer void growth (Raj and Ashby, 1975; Cocks and Ashby, 1982). A void of radius I-/, grows by material flowing along the surface of the void to the point of intersection with the grain boundary and then along the grain boundary, as illustrated in the general picture

Modeling Microstructure Evolution

103

of Figure 1, where it plates out uniformly over an area ~ ( - 1r;),~allowing the applied load to do work. As well as providing material for the growth process, surface diffusion also provides a mechanism for internal rearrangement of material toward a local equilibrium configuration. If 2>, is large, this local rearrangement will occur rapidly and only a slight perturbation from the equilibrium shape, which is sufficient to provide the flux for the growth process, will quickly evolve. In our analysis we can therefore assume an equilibrium profile for the pore, as in the classical model (Raj and Ashby, 1975), and ignore the contribution from the second term in the variational functional of (4.1 1). To simplify the analysis further, we assume that the grain boundary energy is 0, so that the equilibrium shape of the pore is spherical. This is a one-degree-of-freedom problem: The volumetric growth rate of the pore can be directly related to the rate of elongation, Li, of the unit cell of Figure 5 ,

v = 41rr;i,

= xi2ci.

(4.12)

The continuity condition for flux in the boundary is

where r is the radial distance from the center of the pore and j,. is the radial flux. Solving this equation, noting that j,. = 0 at r = I , gives l i ,

j,. = -2r (I-

-

r’).

Substituting this expression for the flux into (4.1 I), noting that

and differentiating the resulting expression with respect to ci to obtain the stationary value, we obtain

where t i , = rill’ is the area fraction of voids on the grain boundary. This expression agrees with the result of Raj and Ashby (1975) as corrected by Raj et al. (1977), although it has been derived in a much more straightforward way. For situations where surface diffusion is rapid, the variational functional only requires a consideration of the flux along the grain boundaries. The above type of

Alan C. E Cocks et al.

104

FIG.6. A closed contour r h around a grain boundary pore used to compute the volumetric growth rate.

analysis can therefore be readily extended to the class of two-dimensional problems examined in Section IV.A, where now the body can contain a distribution of near-equilibrium-shaped pores. Consider the situation where the pores are situated at the triple-points, as is commonly observed in engineering ceramics. We again assume for simplicity that the grain boundary energy is 0, so that the pores adopt a circular profile. The rate of increase of surface area for the pore can then be related to the volumetric growth rate:

There are two contributions to the volumetric growth rate: due to the flux of material away from the pore and as a result of the grain boundary plating process. These contributions can readily be expressed in terms of the degrees of freedom identified above by constructing a closed contour rl, around each pore, as shown in Figure 6 . Then V =

fr

i;nidr+ h

c

j,m,,

boundaries

where n; is the outward normal to the chosen contour, the summation is over all boundaries cutting the contour, and m, is a unit normal along the line of the boundary pointing away from the contour. The velocities u i and fluxes j , can be related to the degrees of freedom [ u ] and thus -

A,&.

= -[U]"F,],

Modeling Microstructure Evolution

105

FIG. 7. ( a ) A macroscopic crack in a polycrystalline material which defomis by Cohle creep. Triple-point voids cxisl along the line of boundaries directly ahead of the crack tip. ( h ) The geometry after an extensive ainount of crack growth.

Alan C. F: Cocks et al.

106 56

-

5452

crack tip position

------

strain rate

-

.-55 50a" .a 48k Y

0

46-

44

-

42-

1 - 1

-/' ____------I

I

where [FYIis a force matrix associated with the change in total surface energy. Substituting this into (4.1 1) and following the derivation of Section IV.A, we obtain

This system of equations can be solved using the same procedures as described in Section 1V.A. Cocks and Searle (1991) have used this method to evaluate the effect of competing length scales on the process of void growth when the spacing of the voids is much greater than the grain size. More recently, Pan and Cocks (1998) have evaluated the process of void growth ahead of a dominant crack. Clearly, the crack tip in this situation does not adopt a near-equilibrium profile (we discuss this more fully in Section 1V.El), but by determining the shape of the crack tip region using the profiles predicted by Chuang et al. (1979) for the same volumetric growth rate they were able to compute the propagation of a crack through an array of triple-point voids lying directly ahead of the crack tip. Computations were performed for regular square arrays of cracks, so that the boundaries shown

Modeling Microstructure Evolution

107

in Figure 7a are axes of symmetry, and for situations where a stress field corresponding to a K I field of prescribed magnitude was imposed on the boundary. A typical profile of the crack and voids after a large extent of crack growth for the former case is shown in Figure 7b and the extent of crack growth as a function of normalized time, as defined in (4.3). is plotted in Figure 8, together with the variation of strain rate for the entire cracked body, normalized with respect to i o , the strain rate of the uncracked body at the same stress. The position of the crack is measured in terms of the number of grain boundary facets that lie along its length. Two processes contribute to crack growth: the diffusional transport of material away from the crack and the growth of triple-point voids. The second contribution results in instantaneous increments of crack growth as the pores link with the crack, indicated by the vertical portions of the curve in Figure 8. It is evident from Figure 7a that toward the end of the life of the component there is significant damage in the uncracked ligament in the plane of the crack, and these voids are all approximately of the same size. Final failure is effectively due to the coalescence of these voids across this minimum section. A full evaluation of these results is provided by Pan and Cocks ( 1998).

2. Fullv Coupled Analysis In this subsection we examine situations in which we cannot assume a nearequilibrium profile for a free surface and the evolving shape is determined by a coupling between surface and grain boundary diffusion. We limit our discussion initially to plane two-dimensional problems, as before. We have already dealt with the boundary diffusion part of the analysis process in Section 1V.A. Here we present an equivalent analysis for surface diffusion. Material rearrangement by surface diffusion results in a change of surface profile. If material is deposited on a section of boundary, it moves in the direction of the outward normal and the continuity equation (3.3) provides a relationship between the flux in the plane of the boundary and the velocity of the surface u,! in the direction of the outward normal n. To model a particular problem, we discretize a section of free surface by a series of straight elements (Figure 9). Having described the initial geometry, care needs to be taken in identifying the degrees of freedom that are used to model the evolution process. Consider the situation represented in Figure 9a in which we associate two velocities with a node at the point of intersection of two elements. In the limit where the normals to the two elements are in the same direction, it is possible for the node to move along the line of the elements without any diffusional rearrangement of material. This is equivalent to a zero-energy mode in a static finite-element analysis and results in a singular surface diffusion viscosity matrix, as defined below. If there is a small angle between the elements, then a much smaller volumetric flux of material is required to translate the node by a

108

Alan C. E Cocks et al.

FIG.9, A section of surface can he discretized by a number ofelements. (a) Two elements with nodal velocities at each node. (h) The degrees of lreedoni and position ofthe Lagrange multipliers for the surface elements.

IWO

given amount in the transverse direction compared with that required to translate the node in the near-normal direction. This results in a low transverse viscosity. This, coupled with the fact that our description of the geometry is approximate (in practice, there are no discontinuities in slope, apart from where a surface is intersected by a grain boundary), results in the transverse velocities dominating the time step and providing only slight modifications to an already idealized situation. It is therefore necessary to constrain the transverse motion of the nodes; i.e., the only degree of freedom associated with a node is in the near-normal direction. The problem now is: How do we define this near-normal direction? There are a number of possibilities, which all produce equivalent results. Pan et al. (1997) choose the velocity to be in a direction normal to a circular arc constructed through the node of interest and its two adjacent nodes. Another approach is to associate it with the direction that bisects the normals of the two adjoining elements. This latter description was used in the simulations described in Section 1V.F. This type of construction is only required for the discretization of a continuous surface. At the point of intersection of a grain boundary with a free surface, there is a distinct discontinuity in slope and two mutually orthogonal velocities are required to model the response. Both types of nodal degrees of freedom are shown in Figure 9b. In practice, the surface and grain boundary tensions at a point where a grain boundary meets the surface must be in equilibrium. This condition is satisfied by

Modeling Microstructure Evolution

109

the variational formalism; i.e., the exact solution to the problem satisfies all the equilibrium requirements. If our initial specified geometry does not satisfy this requirement, then rapid local diffusion will occur to ensure equilibrium of the surface and grain boundary tensions. We need not, therefore, enforce this condition. As with conventional finite-element analyses, when using a discretized description of the problem the optimum solution does not completely satisfy the equilibrium conditions. Similarly, here, the microstructure that evolves may not exactly satisfy the dihedral angle requirement, but by not enforcing this condition we can readily evaluate the appropriateness of our geometric description in the vicinity of a surface/grain boundary junction by comparing the local angle between the intersecting interfaces with that required by equilibrium. We can follow a similar procedure to that described in Section 1V.A for boundary diffusion when examining the contribution of surface diffusion to the variational functional. For a given element, the velocity normal to the boundary at any position can be related to the nodal velocities at its ends. Integration of the continuity condition (3.3) then allows the surface flux to be determined in terms of these velocities and the flux at the center of the element, which represents an additional degree of freedom, giving three in total. Our local flux field does not ensure global conservation of mass and, as in our analysis of grain boundary diffusioncontrolled processes, we introduce a series of Lagrange multipliers, one per node, to enforce the condition that there is zero net flux of material into a node. Also, the rate of change of surface and grain boundary area can be related to the velocities at the nodes. Following the analysis presented for grain boundary diffusion in Section IV.A, we can express the variational principle of (4.1 1) in the following form:

where 3 , ; ' [ K s ] is the surface diffusion viscosity matrix, [Fs] and [ F b ] are the thermodynamic forces determined from the surface and grain boundary tensions, and [C,] is the constraint matrix arising from the conservation of mass condition at the surface nodes. The matrices [ u ] and [A] contain the degrees of freedom and Lagrange multipliers associated with both the surface and the grain boundary diffusion processes. The only common quantities are the Lagrange multipliers associated with the points where the grain boundaries meet a free surface; i.e., the coupling arises from the requirement that the flux into a boundary is equal to the sum of the fluxes away from the adjoining surfaces. These Lagrange multipliers

110

Alan C. I;: Cocks et al.

are effectively the capillarity stresses at these points. All the above matrices are fully defined by Pan et al. (1997). Taking the stationary value of (4.13) gives

This set of equations can be solved using the same standard matrix procedures as for grain boundary diffusion alone and the geometry can be updated by multiplying the velocities by a suitable time step. For large-scale problems a large viscosity/constraint matrix must be inverted at each time step, which can lead to long computer run times. It is evident from the geometric description of Figure 9 that the surface diffusion elements have limited connectivity to other elements, only interacting with two adjoining elements in general. The matrix [&I is therefore a narrow-banded matrix. Under these conditions finite-difference methods can be much more efficient than matrix procedures developed from a variational formulation. Pan and Cocks (1995) describe a combined finite-difference/matrix method for the coupled problem in which the above Euler explicit integration scheme is employed. Recently, Kucherenko and Pan (1998) have described a procedure based on an implicit integration scheme and provided detailed comparisons of the different methods and integration schemes. For problems in which only surface and grain boundary diffusion contribute to the deformation process and in which the grains are assumed to be rigid, these combined schemes are much more efficient than the above matrix method when using the same number of degrees of freedom. When using the numerical schemes developed from the variational principle, very coarse descriptions of the microstructure can be adopted. Even with a small number of degrees of freedom, the computations remain stable, unldce conventional finite-differencemethods. In Section V we illustrate this by considering the cosintering of a row of initially spherical particles by coupled grain boundary and surface diffusion and use the variational principle to model this process in terms of only two degrees of freedom and compare the overall features of the evolution process with that obtained from a finite-difference/matrix method. Numerical procedures based on the variational principle have a distinct advantage over finite-difference methods when other processes or driving forces contribute to the component response. These can produce a stronger interaction between the surface diffusion elements. We consider one such example in Section IV.F, where we take into account the contribution of the elastic deformation of the material to the driving force. It is not immediately evident how finite-difference methods can be adapted to deal with this situation.

Modeling Microstructure Evolution

111

c. I N T E R F A C E REACTIONS In Sections 1V.A and 1V.B we assumed that the surfaces and grain boundaries are perfect sources and sinks for the diffusing species. As noted in Section II.C, some of the work that is available to drive the overall process is required to drive the sources and sinks. In fine-grained materials the kinetics associated with these interfacial processes can control the overall rate of microstructure evolution. The incorporation of interface reactions into the variational principle has been described by Cocks (1992) and Cocks and Pan (1994). No additional geometric information or degrees of freedom are required to model these processes. For a given element, be it a grain boundary or a surface element, we can express the velocity at any point of the element in terms of the velocity degrees of freedom associated with that element. For the purpose of discussion, we concentrate on the situation where interface reactions only operate on grain boundaries; i.e., we consider the surfaces to act as perfect sources and sinks for the diffusing species. The variational functional eq. (3.6) can then be written in the following dimensionless form:

where a! = (iod/vo)'/''' measures the relative importance of the diffusive transport of material and the interface reactions to the rate of dissipation of energy, u, = I U , ~ 1, the magnitude of the rate of plating of material on the grain boundary, and the rate constant uo is related to the constant & defined in (2.8) through the relationship uo = NOR. As before, d is a suitable physical length scale for the problem of interest, which is generally chosen to be the grain diameter. For small a! the kinetics associated with the interface reactions are fast compared to those associated with the diffusional process and only a small proportion of the available work is required to drive the interface reactions. As a result, we can ignore the contribution of the interface reactions to the variational functional and to the overall microstructure evolution process. Conversely, when a! is large the interface reaction kinetics are sluggish, these now determine the overall rate of the process. We can then ignore the contribution of the diffusional flux of ma-

Alan C. E Cocks et al.

112

FIG. 10. A simple rheological model representing the creep behavior of a polycrystalline material. Poc/ and Por are the uniaxial strain rates at a stress 00 in the limits of grain boundary diffusion and interface reaction controlled creep.

terial to the variational functional. Diffusional processes of the type considered here consist of removing material from a section of the grain boundary and transporting it to another part of the body where it is deposited. These are all necessary sequential parts of the same process which must operate at a compatible rate; thus the slowest process determines the overall rate of microstructure evolution. Cocks and Pan (1994) have described this type of response using simple rheological models. For example, the uniaxial strain rate of a polycrystalline material can be represented in terms of the rheological model of Figure 10, where the two parallel dashpots represent the material response in the limits of grain boundary diffusion controlled and interface reaction controlled creep, respectively. The strain rate at intermediate values of a! is given by i. for m = 2, where B = 1 5 . 5 ~ 1= ~ &/&I;, & is the uniaxial strain rate in the limit of grain boundary diffusion controlled (Coble) creep, and Boi is the uniaxial strain rate in the limit of interface reaction controlled creep.

Modeling Microstructure Evolution

113

Following the analyses of Sections 1V.A and IV.B, we can express the variational functional in terms of the velocity and flux degrees of freedom and the nodal constraints, represented by the Lagrange multipliers. The values of these quantities, which provide a stationary value of the functional, satisfy the following set of nonlinear simultaneous equations:

where 1K , 1 is the interface reaction viscosity matrix, which is a function of the degrees of freedom [ u ] . Cocks (1992) describes a Newton-Raphson method for solving this system of equations. Cocks and Pan ( I 994) have used this method to evaluate the process of void growth in fine-grained materials and to determine the nature of the stress and velocity fields that develop ahead of a crack tip in a creeping material, and how these depend on the parameter a. Consider the situation where a body is subjected to a single force F and u is the velocity in the direction of F . For the purpose of discussion, we only consider the dissipative processes of grain boundary diffusion and interface reactions and ignore all other contributions to \I, and G. In the limit of grain boundary diffusion controlled microstructure evolution, i.e., when a is small, the variational functional of (3.6)becomes

fl,.=

L,,5

1;

7

ju ja d

A -FC.

This functional can be minimized by minimizing the distance that the material diffuses, subject to the constraint that the resulting velocity field is compatible. In polycrystalline materials, a compatible mechanism can be obtained by the material diffusing a distance of the order of the grain diameter. In all the problems that have been analyzed, we observe that material never diffuses distances greater than this, even when the body contains large defects. For example, in the crack problem of Figure 7 the high plating rates in the vicinity of the crack tip are accommodated by Coble creep of the surrounding cage of grains and material that flows into this region from the growing crack is plated onto the first boundary facet. In the limit in which the interface reactions control the rate of microstructure evolution, i.e., for large a,the variational functional of (3.6)becomes

For the situation where the body contains defects, in the form of either voids or cracks, the functional is now minimized by material flowing away from the defects and plating uniformly over all available facets that allow the applied load

114

Alan C. E Cocks et al.

to do work. For the crack problem of Figure 7, we now find that material from the growing voids and crack plates out uniformly on all the grain boundary facets which make an angle of 30" with respect to the line of the crack; i.e., material flows from the surface of the crack to the extremes of the repeating cell. For a macroscopic crack this limit will never be achieved in practice. In fine-grained materials the velocity and flux patterns are determined by both dissipative processes. Cocks and Pan (1994) demonstrate that the distance that material diffuses depends on the value of a. They demonstrate that, on average, material diffuses a distance that scales with ad if a > 0.5, while the diffusion distance is of the order of d / 2 for smaller values of a. Thus, if a fine-grained material contains a macroscopic crack extending over hundreds of grains, a diffusion zone of mean dimension ad forms ahead of the crack tip. Material that diffuses away from the growing crack remains within this zone and the resulting deformation is accommodated by interface reaction controlled creep of the surrounding network of grains.

D. COUPLED G R A I NBOUNDARY D~FFUSION A N D SELF-DIFFUSION In the previous subsection we considered the situation where two sequential processes contribute to the material response. Here we consider a simple situation where two parallel/competing processes occur. We simply examine the form of the variational principle and do not attempt to develop numerical procedures from the variational formalism. Consider, for simplicity, the situation where only lattice and grain boundary diffusion contribute to the component response, so that, by default, there are no free surfaces in the body. The variational principle then becomes, in dimensionless form,

(4.17)

where L' $ = Vld/'D/,. The lattice and grain boundary fluxes are, in general, independent of each other and combine to produce the displacement rate l i i , particularly at extreme values of 81.In situations where V / is of the order of unity, flux

Modeling Microstructure Evolution

115

patterns can develop where the lattice and grain boundary flux fields are incompatible when taken in isolation. The actual rate of energy dissipation is then faster than that determined by simply summing the contributions from the two mechanisms operating in isolation. This can be readily illustrated by considering the situation where the applied surface tractions T, are specified. The rate of energy dissipation associated with the exact field is then

In the limit where grain boundary diffusion dominates,

where the superscript b indicates rates determined by only considering grain boundary diffusion. Only considering material rearrangement due to lattice diffusion. we find

where the superscript I indicates the fields determined in this limit. According to the variational principle, the value of I?,. determined from the full coupled analysis is less than or equal to the value obtained by adding the two extreme mechanisms. Mahng use of the above relationships, we find

If the body is subjected to a single load on its boundary and the rate of change of grain boundary energy is small compared to the external work rate, then the displacement rate in the direction of the applied load is greater than the sum of the displacement rates for each mechanism operating in isolation. For small 31it is evident that the minimum value of fi,. occurs when most of the displacement rate is due to grain boundary diffusion, so that the second term of the functional is comparable to the first. Similarly, if 2)/is large, lattice diffusion dominates the response. This is an obvious example, but it illustrates the way in which the variational principle can be used to determine the relative importance of different mechanisms in a given situation. We consider other examples of this in subsequent sections. Numerical and approximate procedures for situations where lattice diffusion is the dominant mechanism have been described by Cocks (1996).

116

Alan C. F: Cocks et al.

E. G R A I NB O U N D A RMIGRATION Y In this section we examine how the variational principle of Section 111 can be used to model situations in which grain boundary migration contributes to the way in which a microstructure evolves. We start by considering the process of normal grain growth in Section IV.E.1, in which grain boundaries move under the action of grain boundary tension. In the process of normal grain growth described here, the only factors that determine whether a grain grows or shrinks are the local topology of the grain network and the grain boundary geometry, which wants to satisfy surface tension equilibrium. Abnormal (or secondary) grain growth, which occurs when some grains have a preferential growth rate due to some energetic advantage over the others, is considered in Section IV.E.2. In Section IV.E.3 we examine situations in which grain boundary migration combines with other diffusional processes to determine the way in which the microstructure evolves.

1. Normal Grain GI-owrh We consider the case here of two-dimensional normal grain growth (Cocks and Gill, 1996; Gill and Cocks, 1996; Du et al., 1998). This removes some of the geometric complexity of working in three dimensions and the results have the benefit of being easier to visualize and interpret. Pseudo two-dimensional grain structures do exist in thin films in which the mean grain diameter of the grains is larger than the film thickness. Such grain structures are said to be columnar and can be represented by a two-dimensional plan view of the microstructure. To represent such a microstructure mathematically, we need a geometric description for the grain boundaries and a topological description of how these boundaries are organized to form the grain network. In this subsection we assume that the only dissipative mechanism is grain boundary migration. Then, from (2.13), for a thin film of thickness h ,

where Lr, = A/,/ h is the total grain boundary length, and the driving force for microstructural evolution is the reduction in the total grain boundary energy

6 =h

IL

Y/,K U,,dS,

h

where K is the curvature of a boundary in the plane of the body. In addition, equilibrium of the grain boundary tensions must be considered. Assuming that the grain boundaries are chosen to be smooth and continuous, the

Modeling Microstructure Evolution

117

FIG. I I , For uniform grain boundary energy per unit area y/,. ( a ) an equilibrated grain boundary junction and (b) a grain boundary junction in which the grain boundary tensions are not in equilibrium leading to a net force F.

only location at whch equilibrium must be enforced is where the grain boundaries meet. In two dimensions, grain boundary junctions of more than three boundaries are assumed to be unstable and, consequently, grain boundaries are said to meet at a triple-point. To satisfy equilibrium, the angle between any two boundaries meeting at a triple-point must by 120" as shown in Figure I la. The microstructural description can be chosen so that this condition is automatically satisfied. Alternatively, this condition can be imposed in a weak form through the variational formulation. Consider the situation shown in Figure 1l b of a triple-point that is not necessarily in equilibrium. There is a net force acting at the triple-point and, therefore, for uniform grain boundary energy yt,,

i=l

where si are the unit tangents to the grain boundaries at the triple-point. If the triple-point is in equilibrium, then F = 0. If the triple-point has a velocity v, then the rate of work done is F v and hence the variational functional can be written in dimensionless form as

-

all triple points

(4.18) potnis

Alan C. l? Cocks et al.

118

FIG. 12. Thc cubic polynomial grain boundary description. The boundary is completely defined by the position ( - 4 - 1 , y1, .r2, y2) and orientation ( H I , 8 2 ) of its two end points. assuming a profile of the forin uI(.y)= rrs' /xs2 C S t / .

+

+ +

where n, = n,/cq)uohRo, F = F/oohRg, V,, = v,,/uo, S = s/Ro, and k = K R O with , 00 = yh/Ro, u g = Ml,oo, and Ro representing the initial mean grain size. We describe a given polycrystalline microstructure by the position of its triplepoints and the orientation of the grain boundaries at these triple-points. Hence we wish to choose a grain boundary description using a smooth function that is completely defined by the position and orientation of its end points but also one that allows their position and orientation to vary independently of one another. The simplest function of this type is a cubic polynomial. A typical boundary is shown in Figure 12. The rate of evolution of this boundary is completely described by the four translational and two rotational velocities of its two end points, which combine to provide the elemental contribution to the degree-of-freedom matrix [ u ] . By considering an infinitesimal temporal displacement of this boundary subject to these end point velocities, it is possible to find a linear relationship between the normal velocity of the boundary u,, at any point and the end point velocities, i.e., u,, = [B,,,l[u],where [B,,] relates the velocity at a material point to the nodal velocities. The variational functional can then be discretized in terms of the triplepoint velocities given the current geometry of the microstructure. The consequent minimization of the variational functional therefore gives the actual triple-point velocities in a similar way to that described in Sections 1V.A and 1V.B for situations in which grain boundary or surface diffusion are the dominant mechanisms. This allows the microstructure to be updated over a suitable time increment. The construction of the initial grain structure is described by Gill and Cocks (1996). They model a large periodic array by identifying a repeating cell and connnecting the boundaries on one side of the cell to those on the opposite side. However, this periodicity means that the viscosity matrix is sparse (cannot be diagonalized) and hence it is computationally more demanding.

Evolution ofthe microstructure Inspection of (4.18) shows that the last term, introduced to ensure equilibrium at the triple-points, acts as a driving force for

Modeling Microstructure Evolution

(a)

(b)

119

(c)

F I G . 13. Dcgrees offreedom o i a triple-point for grain growth models: (a) Case 1, (b) Case 2, and (c) Case 3.

grain growth even when all the boundaries are straight (i.e., G = 0). Therefore there are two terms that contribute to the driving force. This suggests three different models by which the microstructure could be allowed to evolve:

Case 1 : Straight boundaries and hence not in equilibrium at triple-points ( G = 0 but F # 0). Case 2: Cubic boundaries fixed to be in equilibrium at triple-points (G # 0 but F = 0). Case 3: Cubic boundaries not geometrically constrained to be in equilibrium at triple-points (G # 0 and F # 0). As shown in Figure 13, the first case is the most restrictive model as each triplepoint only has two degrees of freedom (both translational), followed by Case 2 which has three (one rotational and two translational), and, finally, Case 3 which has five (three rotational and two translational). As grains are growing and shrinking, it is necessary to incorporate the topological changes that allow this to happen into the updating procedure. These topological processes can be modeled by the combination of the two types of events shown in Figure 14 (Ashby and Verrall, 1973). The first of these, the grain boundary switching event, occurs when a boundary connecting two grains A and B becomes very small and the grains A and B come into contact. This is modeled by reorienting a boundary when it shrinks below a critical length. The second, three-sided cell removal, is simply the contraction of a three-sided cell to a single point at a grain boundary junction. This is modeled by removing three-sided cells when their area falls below a critical size. Combining these processes allows grains with any number of sides to be removed. In the numerical simulations it was assumed that the grain boundary energy y/, is constant. The solution is found to be independent of the film thickness h and the time is normalized with respect to the only material parameter MhaoRo =

120

Alan C. E Cocks et al.

(b)

FIG. 14. Schematic diagram of the two types of topological changes: (a) neighbor switching and (b) disappearance of a three-sided grain.

Snapshots of a typical Case 2 simulation starting with 1024 grains are shown in Figure 15 (Gill and Cocks, 1996). Simulations of the simpler Case 1 model conducted by Du et al. (1998) possess the same general features of the Case 2 model. For a given grain structure the Case 1 model requires fewer degrees of freedom to describe the microstructure, resulting in a reduction in the total time required for a simulation. In each case, however, as in all the numerical simulations considered in this paper, a system of N, (where Nt is the number of degrees of freedom used to describe the system) linear simultaneous equations is solved at each time step. Du et al. (1998) have used the variational principle to develop an approximate nodal model for an array of straight-sided grains. The force matrix, with each component representing the force on the associated triplepoint, is determined as for the Case 1 model. A viscosity matrix is determined for each node by arbitrarily associating half of each grain boundary facet with its two end nodes and assuming that the normal velocity along the length associated with a given node is constant and equal to that at the node. The velocity of a given node is then determined by solving two local linear simultaneous equations; thus the rate of change of microstructure is determined by solving Nt / 2 pairs of equations. Although approximate, this procedure is much quicker than inverting the full viscosity matrix at each time step. Simulations using this method differ in detail from those produced using the full variational method, but when the results are examined in terms of various macroscopic parameters there is good agreement between the two sets of simulated results. YbMh.

Modeling Microstructure Evolution

T i m =OOO

Numer of cell 4 0 2 4

T i m =I95

Number of cel?, -707

Time =a54

121

Nmber of cdk =269

f)

I

Time =4 21

Nmber of cels -464

Flc;. I S . Thc grain growth evolution of a network of p i n s as a function or noimalized time T = t/y/,M/,:(a) starting with 1024 grains, (b)-(rS some grains grow at the expense of others causing sonic of ihern to disappeai-. Reprinted front Gill and C o c k (1996) with permission of Elsevier Science.

When a similar simulation to that shown in Figure 15 was conducted for the more general Case 3 (Gill and Cocks, 1997b), it was found that the grain growth occurred roughly 5% quicker but the general characteristics of the evolution were exactly the same. It is thought that the growth was slightly faster as the additional

122

Alan C. E Cocks et al.

degrees of freedom allow a slightly more optimal energy reduction path to be followed. However, the principal result from this study was the fact that the computational time required was 10 times that for the simpler Case 2 study as the more complex model proved to be less stable. This clearly demonstrates that one has to be careful in choosing the degrees of freedom for a given problem. If the problem description is not given enough degrees of freedom to encompass the physical behavior of the system, the results will be misleading. However, if the problem is given a lot of freedom, it may become less stable, i.e., the functional minimum becomes shallow, and it may not become possible to obtain a solution at all. To get the full value out of these simulations, it is necessary to investigate the evolution in more detail and extract some simple laws or parameters that more readily encapsulate the characteristics of the simulation. There are many statistical methods and theories to monitor and predict the kinetic and topological evolution of cellular networks (Atkinson, 1988; Gill and Cocks, 1996). The most widely used model for grain growth kinetics is the parabolic grain growth equation proposed by Burke (1 949): R” - R; = k t ,

(4.19)

where n takes the value of 2 and is known as the grain growth exponent, R is the mean grain radius at time t , and k is a rate constant. As the mean grain area (taken to be R 2 ) is inversely proportional to the number of grains N , the grain growth exponent obtained from the Case 2 simulations can be readily calculated. In this case it was found to be n = 1.92 in reasonable agreement with (4.19). In real material systems, however, n is found to be between 2 and 4. This discrepancy will be discussed further in the section on abnormal grain growth. Rhines and Craig (1974) introduced the concept of a sweep constant in an attempt to bridge the gap between the topological and kinetic models of grain growth. Their definition was later revised by Doherty (1975) to define the sweep constant O* as “the number of grains lost when the grain boundaries sweep through a volume equal to the mean grain volume.” The numerical results were found to justify Doherty’s definition with a constant value of O* % 1.32. This definition is important for the development of a simplified variational model. Simplified variational approach The numerical results demonstrated that there is significant evidence that a number of parameters exist that can accurately represent the characteristics of normal grain growth evolution. These are fully described by Gill and Cocks (1996), who demonstrate that it is possible to represent the observed self-similar evolution of the system by a single variable. One choice for this degree of freedom is N ( r ) , the number of grains in the network as a function of time r .

Modeling Microstructure Evolution

123

Let the constant Ao be the total cross-sectional area of the unit cell and let A be the average area swept per grain boundary per unit time. Then the number of grain boundaries in the unit cell is 3 N (as the average number of sides in a trivalent cellular network is necessarily 6), the number of grains lost per unit time is -N, and the mean grain size at time t is A o / N . Thus the sweep constant is given by

If the average grain boundary length is L and its average absolute normal drift velocity is u , , ~ ?then . A = u,,,,? L . The total grain boundary length L, = 3 N L is empirically found to be Lt = c2N1/’, where c i = 2 f i a ’ A o and (Y = 0.96. Combining these expressions yields

A

utlln

N

=---cI- N I / Z ’

L

where c: = A 0 / 2 f i ( a e * ) ’ . The variational functional (4.18) can be expressed in terms of these average quantities as

where u,~,.is the root mean square of the normal drift velocity, which has been written as u:,. = q(u,,,,,)’ and q > 1 is an unknown constant. This functional is a minimum when d & / d ~= 0, which gives, upon rearrangement and integration,

,/m,

As the mean grain size R = it can be seen that the above is similar in form to (4.19) with k = f i ( ( ~ O * ) ’ y [ , M b / q .Hence we have shown that the variational principle predicts a parabolic grain growth law without making any assumptions about the grain boundary dynamics except that 8* is a constant. We have also derived an expression for the time constant k in terms of 8* and an unknown parameter q. From the numerical results we find that q = 2.76 and 8* = 1.32, which predicts that k = 1.01, in good agreement with the actual value of k = 0.98. The important point to note about the above analysis is that as well as providing the theoretical background for the development of the numerical procedures the variational principle also allows the results of the simulations to be evaluated and

Alan C. F Cocks et al.

I24

a simple macroscopic law to be developed. We extend this concept further in the following subsection. 2. Abnormal Gmin Growth The evolution of polycrystalline thin films can be strongly influenced by the substrate on which the film is deposited. In general, the interfacial energy between the film and the substrate is lower if a grain has a similar crystallographic orientation to the substrate than if it does not. This variation in the interfacial energy between grains means that the energetic conditions for the growth of some grains are more favorable than those for the growth of others. This leads to abnormal (or secondary) grain growth (Thompson, 1990). In this case, unlike in normal grain growth, no time-invariant grain size distribution exists. Abnormal grain growth is characterized by the development of a bimodal grain size distribution in which abnormal grains (those with a low film-substrate interfacial energy) generally grow at the expense of so-called normal grains (those not preferentially favored for enhanced growth). The variational approach is flexible enough to incorporate many driving forces and dissipative mechanisms. The functional described here is similar to that for normal grain growth (4.18) but with an additional driving force. For the purposes of clarity and illustration, it is assumed here that the grain boundary energy per unit area, p,,is constant and that there are only two different interfacial energies. The difference between these energies, Ax,,, is the driving force for abnormal grain growth. The variational functional is therefore (Gill and Cocks, 1998)

h

u i ds

+ hy17

L,,

+

K U ds ~ ~

-

l),

F v f Ayln grain boundary junctions

u,, ds

or, in dimensionless form, using the same normalizations as for normal grain growth, (4.20) all triple points

where the final term is the rate of change of Gibbs free energy due to the change in the total energy of the film-substrate interface and the parameter (4.21)

Modeling Microstructure Evolution

125

FIG. 16. The coalescence of two abnormal (light) grains in a matrix of normal (dark) grains.

where Ro is a reference length taken to be the initial mean grain radius. This parameter has been proposed elsewhere (Floro and Thompson, 1993) and represents the propensity of a system toward abnormality.

Evolution of the microstructure The initial microstructure is generated as for normal grain growth. A certain percentage of the grains are randomly chosen to be abnormal, i.e., to have a film-substrate energy that is A x n less than the others. The boundaries between abnormal grains are assumed to have a very low energy as they have a similar crystallographic orientation. Consequently, such boundaries are simply removed from the simulation. The process of two abnormal grains coming into contact is therefore one of coalescence, as shown in Figure 16. Once a boundary is removed, the two boundaries that remain at a junction are no longer in equilibrium. The microstructure is therefore allowed to restore itself to equilibrium under the action of the grain boundary tension by assigning such boundaries independent rotational freedom. Otherwise, the triple-points are chosen to be of the type Case 2 (see Section IV.E.l) to facilitate the computational process. Typical values for the relevant material parameters are approximately A x n = 0.05 J/m2, yb = 1 J/m', h = 50-500 A, and Ro = 500 A (Floro and Thompson, 1993). This gives a range for realistic values of 20of 0 5 Zo 5 0.5, where

126

Alan C. E Cocks et al.

Modeling Microstructure Evolution

127

128

Alan C. E Cocks et al.

Modeling Microstructure Evolutioii

129

20 = 0 corresponds to normal grain growth. The results of two typical simulations for the same initial microstructure and ZO = 0.047 and Zo = 0.47 are illustrated in Figures 17 and 18, respectively. Gill and Cocks (1998) have fully evaluated the results of these simulations. Here we examine some of their major observations. It is immediately evident from Figures 17 and 18 that for small 20 there is significant normal grain growth before the growth of the abnormal grains takes over and consumes the entire film, while for the larger value of ZO there is limited normal grain growth and that preferential growth of the favorably orientated grains occurs early in the evolution process, with the final size of the normal grains only slightly larger than at the start of the simulation. Also, when ZO is small a significant number of the preferentially orientated grains are lost during the early stages of grain growth, while the majority of these grains remain when Zo is large and contribute to the abnormal grain growth process. After a time t, the abnormal grains consume all the normal grains to produce a single crystal. Over the range of conditions considered in the simulations, this time is well represented by the following empirical relationship:

t- =

.

g ~

(Zo+C.)’

where c and g are dimensionless material constants, which depend on the initial abnormal grain area fraction f . For ,f = 4%, c = 6.3 x lo-’ and g = 6.4. The area fraction of abnormal grains at a given time, , f ( t ) , is of practical interest to experimentalists. Figure 19 illustrates the temporal evolution of this quantity for the two numerical simulations mentioned above. As we shall see later, it is of great interest to know the length of the grain boundary interface between the abnormal and normal grains, L(,. Von Siclen (1996) predicts that L,, = ( 1 - f’)L,‘,X, where L:,X is called the “extended” length of the abnormal-normal interface and represents the length of the boundary in the absence of grain impingement. Using the Kolmogorov-Johnson-Avrarn-Mehl equation (see Von Siclen, 1996), f = 1 - exp(-fex), where f“ is the extended abnormal grain area fraction, and, assuming that grains can be represented as circular in the absence of impingement, we find L,, = 2(1 - flJ-arN,,Aolog(l

-

f),

(4.22)

where the total area under consideration, Ao, the initial number of abnormal grains, N,, and a , the fraction of abnormal grains that continue to grow and do not shrink, are all constants. The parameter a can be estimated from Hillert’s ( 1965)

130

Alan C. E Cocks et al.

F I G . 18. Snapshots of the abnormal grain growth evolution of a polycrystalline thin film microstructure for Zo = 0.47 at normalized time: (a) t = 0.0 s, (b) t = 1.4 s. (c) r = 2.8 s, (d) r = 4.8 s. (e)r = 6.7 s. and (f) I = 8.4 s.

Modeling Microstructure Evolution

131

a

Modeling Microstructure Evolution

133

grain growth law with an additional driving force term due to the interfacial energy difference. A grain of radius R, is assumed to grow initially if (4.23)

+

i.e., when R,/Ro > 1/(1 2Zo). Assuming the initial grain size distribution can be represented by the steady-state grain size distribution for normal grain growth proposed by Hillert (1965), we can find the fraction of grains that satisfy this condition 3

a!

+4zo

= (iTG)exp(&).

This predicts that a! = 0.50 for ZO = 0.047 and a! = 0.85 for ZO = 0.47. For the simulation shown in Figure 17 for Zo = 0.047, it was found that only 30% of the abnormal grains survived (given the initial random microstructure shown in Figure 17a) so a value of a! = 0.30 is used in subsequent calculations. The model is much better for larger ZO = 0.47, when a much larger fraction of the favorably orientated grains survive and contribute to the abnormal grain growth process. Simplijied variational approach To develop a macroscopic model of abnormal grain growth, we return to the variational principle of (4.20) and, following the analysis of normal grain growth, express the material response in terms of a limited number of state variables. There is one driving force for normal grain growth and the state of the microstructure at any point can be represented by one variable, N i l ,the number of normal grains. For abnormal grain growth there are two driving forces and so an additional variable is required. We choose the abnormal grain area fraction f . Gill and Cocks (1998) have expressed the variational functional of (4.20) in terms of these two state variables and by determining the values of f and N,,that minimize the functional obtained evolution laws for both J' and N i l . Here we concentrate on the evolution of the area fraction of normal grains: (4.24) where the parameter K",, is the increased curvature of the abnormal grain boundaries over that of the normal grain boundaries due to changes in the grain boundary geometries from geometrical constraints. This is the net mean curvature of the abnormal grain boundaries, i.e., a directional mean curvature taking the curvature to be positive if its center lies on the normal grain side of the abnormal-normal grain interface. It has been found from the numerical results that K,,, quickly attains a reasonably constant value, once any abnormal grains that are not going to

Alan C. E Cocks et al.

134 0.9

0.8 0.7 f 0.6 0.5 0.k

-ZO = 0 047 - - - Z0=0.47

0.3

01

0

I

0.1

I

0.2

I

0.3

I

0.4

I

0.5

I

I

0.7

0.6

I

I

0.9

0.8

F I G .i 9. The abnormal area fraction f as a functioii of time I normalized with respect to the total evolution time r for 20 = 0.047 and Zo = 0.47.

survive have disappeared, until the very final stages of growth, during which the net curvature can increase. Consequently, using (4.22), assuming a constant value of K(!,, , and integrating (4.24) after separation of variables, we obtain f = 1 - exp(-h2(t - t i ) 2 ) ,

(4.25)

where h=

M O Y O ~ ( z f0 ROKori)

4%

,

tl

=

J-

- fo)

h

and fo is the abnormal grain area fraction at time t = 0. Equation (4.25) is the classical expression commonly used to obtain a fit to experimental data (Thompson, 1985).However, the asymptotic behavior of (4.25) as f -+ 1 is not correct, i.e., f + 1 only as t + 00, although the majority of the evolution is faithfully represented. The quality of the fit is determined by the accuracy of the model for the abnormal grain boundary length, eq. (4.22). This equation provides a good representation for 20 = 0.47 and hence eq. (4.25) provides a good fit to the numerical results of Figure 19 except for the final stages of abnormal grain growth, i.e., f > 0.9. The fit is not as good for 20 = 0.047 when there is a much stronger interaction between the normal and abnormal grain growth processes.

Modeling Microstructure Evolution

135

We have investigated here the case of a thin-film microstructure in which there are two different film-substrate interfacial energies. A more general formulation for multiple interfacial energy systems could be obtained using a variational approach with multiple state variables. This is discussed more fully by Gill and Cocks (1998).

3 . Coupled Grain Boundap Migration and Surjiace Diflusion In the above analysis we considered a two-dimensional plan view of the surface of the thin film when modeling the process of abnormal grain growth. In practice, the surface profile is not flat. Grooves can develop where the grain boundaries meet the free surface, with the depth of these grooves and the surface profile depending on the velocity of the migrating grain boundary. A two-dimensional analysis of this situation has been presented by Suo (1996) and Mullins (1958). Recently, Pan et al. (1997) modeled this problem by combining the surface diffusion relationships of Section 1V.B with the variational functional of Section IV.E.2. Using the same normalizations as in Section IV.E.2, the variational functional becomes

poiII t s

where 8f,,,s = D , s , / M ~ Rand i Ps = yr/y/,, which we assume to be constant for the entire surface. If fi,,,,, is small, then the surface diffusion hnetics are sluggish compared with the rate of grain boundary migration. The surface therefore remains flat as the boundaries migrate and the analysis of Section IV.E.2 is valid for the full three-dimensional situation. If Bfii,, is large, significant grooving of the grain boundaries will occur before the boundaries have had the chance to move, which will effectively pin the boundaries, preventing their migration. When is of the order of unity, both surface diffusion and grain boundary migration determine the way in which the microstructure evolves. Pan et al. (1997) have used the variational principle to analyze the twodimensional situation considered by Mullins (1958) and Suo (1996), which is shown in Figure 20a, and consists of a single migrating grain boundary. The only length scale in this problem is the film thickness h , which we take as our normalizing length, rather than the grain size Ro. The steady-state profile is shown in Figure 20b for fi,,,,= 0.096, 20 = 0.5, and y, = 1.0 and compared with the predictions of Suo and Mullins. These values were chosen to violate the small-

136

Alan C. R Cocks et al.

---_-Finite element solution. Steady state solution by Mullins and Suo.

(b)

FIG.20. Thermal grooving at a migrating grain boundary. (a) The initial geometry. (b) The steadystate solution compared with lhe analytical solution of Mullins (1958) and Suo (1996). Reproduced from Pan e/ a/. (1997) with permission of The Royal Society.

slope assumption of the analytical solution (in the analysis of Suo and Mullins, it was assumed that the curvature of the surface is approximately given by the second derivative of height with respect to distance along the film). Despite this, there is good agreement between the two sets of results. A fuller evaluation of the predictive capability of this variational principle is given by Pan et al. (1997). Full three-dimensional calculations of the evolution of the grain structure in a thin film are yet to be undertaken.

F. EFFECTOF CHANGES IN ELASTICSTORED ENERGY ON MICROSTRUCTURE EVOLUTION In this section we add the contribution of changes of elastic stored energy to the variational principle. Two methods of analysis are described, which make use of boundary element and finite-element methods to determine how the Gibbs free energy changes as the microstructure evolves. We use the first of these approaches to evaluate the diffusive growth of a grain boundary crack. This is presented in Section 1V.F. 1. The second method is described in Section IV.F.2 and is employed to evaluate the growth of instabilities in thin films. 1. Difusive Growth of an Intergranular Crack in an Elastic Bicrystal

In Section 1V.B we examined the classical model of Hull-Rimer void growth in which it is assumed that surface diffusion is sufficiently rapid for the pores to

Modeling Microstructure Evolution

137

F I ~ ;2. I , Hall geometry of the intergranulai- crach problem for a linite crack of halt' width ( I . The stress licld due to the prewnce of thc crack o<.(.v)is singular at the crack tip. Material elastically wcdgcd into the grain boundary causes ;in opening displacement 6( Y ) and alters the stress distribution along the boundary. Thc displacemeni utlopta a confguratioti such that the stress singulaiity at thc crack tip disappears and the resulting stress normal to the houndary O , ~ ( . I) is finite.

maintain an equilibrium profile as they grow. In this section we examine the situation where both surface and grain boundary diffusion processes determine the rate of void growth (Chuang ef ul., 1979; Thouless et al., 1983).In a rigid material, any matter removed from the crack must be uniformly plated over the grain boundary. This situation has been investigated by Pharr and Nix (1979). In an elastic material, however, this constraint does not apply as nonuniform plating of material over the grain boundary can be accommodated by elastic defoimation of the matrix. The elastic stored energy associated with this wedging of material into the boundary depends on the current plating distribution (or opening displacement) over the entire length of the boundary, and, consequently, it depends on the entire history of plating on the grain interface. Chuang ( 1982) analyzed the steady-state propagation of a semi-infinite crack in an elastic bicrystal. This necessarily assumes that a steady-state plating distribution exists and is reached within a finite time. The method presented here is a numerical model for the growth of a finite crack which is dependent on the loading history, and hence the plating history, of the specimen. Consider an intergranular cavity of finite length 2u at the center of a grain boundary of length 2w,as shown i n Figure 21. The height of the cavity is assumed to be small with respect to its length, i.e., it is a cracklike cavity, and plane strain conditions are assumed. The x axis is taken to coincide with the grain boundary with the origin at the center of the crack so that the ends of the crack are at

138

Alan C. I? Cocks et al.

x = f a . The remotely applied stress

acts normal to the grain boundary and is assumed to be constant. The dissipation rate in the variational functional contains terms due to grain boundary and surface fluxes 0 ,

\Ir=\k,+\Irb=

and is accompanied by suitable Lagrangian terms to ensure matter conservation as described in Sections 1V.A and 1V.B. The driving force arises from the rate of change of the surface area of the grain boundary A b and crack surface A, and the rate of change of total potential energy of the system A so that

G = y S A &+ Y,,Ah

+ A.

The variation of elastic stored energy in diffusive cavity growth has been considered by Chuang and h c e (1981). Assuming that the strain energy density generated is low, even near the crack tip, energy variations due to the flux of strained material can be considered negligible. The rate of change of potential energy is then simply the difference between the rate of change of elastic strain energy due to removal of material to extend the crack and work done rearranging material on the grain boundary. If the stress normal to the boundary at any instant is a,,(x) and the opening displacement of the boundary is 6(x), then Chuang and Rice (1981) demonstrated that

(4.26) Chuang (1982) proposed that the boundary opening displacement could be represented by a continuous array of infinitesimal dislocations. This approach will be adopted here and is outlined below for the case of a crack of finite length. Consider an edge dislocation with Burgers vector b a distance u in front of the crack. The misfit stress at a position x on the grain boundary associated with this dislocation is (Luo, 1978)

a;n ( u , x-> =

bE

4 ~ ( -l U * ) ( . X - U )

for a material of Young’s modulus E and Poisson’s ratio u. The formation of a dislocation is due to a change in the opening displacement of the boundary and, following Chuang (1982), we write b = -d8/dx = -S’(x). Consequently, the stress normal to the grain boundary is the stress due to the presence of the crack

Modeling Microstructure Evolution

139

F I G . 2 2 . Discretization of the crack problem. Thc evolution of the geometry is described by the nornial velocities to the crack surface 1 1 ; and the opening displacement ratcs 8, at given points on the grain boundary.

a,.(x) and the sum of the stresses due to the array of dislocations

where a,, is the mean stress acting normal to the boundary. The stress distribution can be determined at any instant given a knowledge of the opening displacement distribution 6(x). Therefore the current state of this system can be represented by a discretization of the crack surface profile and the grain boundary opening displacement, as shown in Figure 22. Modeling the motion of a free surface by surface diffusion has been discussed previously in Section 1V.B. The opening displacement distribution is described by a number of linear elements which are defined by their end points 6,. However, the element at the crack tip requires further consideration. The crack tip opening displacement adopts a configuration so that the stress field associated with it interacts with the stress field due to the presence of the crack to remove the singularity in the stress field at x = a. This suggests that the stress field due to the opening displacement in the locality of the crack tip

Alan C. I? Cocks et al.

140

must also be singular. Inverting the Cauchy integral in (4.27) gives

Assuming that the stress at the crack tip is finite, one can evaluate this integral for x in the locality of the crack tip to give S’(X)

=

I+

4(1 - u2) X - U ITE all( a )In I T

finite terms

for (x - u ) / a << 1. Thus the opening displacement at the crack tip is described by

where a is an unknown constant dependent on the stress at the crack tip, a,( a ) , and the opening displacement at the crack tip, 60, is chosen to ensure continuity between elements. The stress distribution along the grain boundary can now be evaluated in terms of the discretized opening displacements 6; from (4.27). This yields an expression of the form

where A and B are constants and fi,, (1)is singular to the order ln(x). When evaluating the rate of change of elastic stored energy (4.26), the integral of the firstorder singular stress component O;,(x) is finite. The first term represents the stress singularity at the crack tip and can be removed by prescribing a = - B / A . This determines the value of stress at the crack tip. The constant B depends on an integral over the entire opening displacement distribution and A depends on the size of the opening displacement tip element. It is interesting to note that the piecewise linear opening displacement distribution yields a finite elastic stored energy, but not a finite stress. Also, a piecewise step function description for the opening displacement, such as that utilized by Thouless et al. (1983), does not yield a finite stress or finite elastic stored energy. The stress presented in the numerical results is the average integral of the stress over a single element. The variational functional can be formulated in terms of the degrees of freedom of this system, the normal velocities to the crack surface u; , and the opening displacement rates & at the discretization nodes, with appropriate Lagrangian terms to ensure mass conservation between the elements; see Figure 22. The geometry can then be updated using a suitable time step. Implementation and validation

Modeling Microstructure Evolution

141

of the method is presented by Gill et al. (1998). Here we present the results of some simulations which allow a direct comparison to be made with the results of Chuang (1982). Rather than write the variational principle in dimensionless form, it proves more convenient here to retain the full expression for the variational principle and employ material parameters for a representative high-temperature ceramic. Following Chuang ( 1982), we use the following material parameters, which represent the response of Salon at 1400°C: y5 = 0.75 J/m’, yl, = 0.375 J/m2, E = 200 GPa, w = 0.27, D, = 2.52 x m4kg-’s, and DI,= 2.75 x m4kg-’s. The initial crack was propagated from an equilibrium cavity of half length a0 = 1 p m on a grain boundary of total half length w = 100 pm. This is not a cracklike cavity, but a cracklike profile develops almost immediately. The remote applied stress was held constant at a, = 100 MPa and so the mean boundary stress a,, increases with the crack length. Physical quantities are normalized with respect to a reference length ao, the initial radius of the void, a reference stress ,,a and a reference time to = ai/Dha,. In this case ro = 420 days. This is large because the applied stress is quite low. Chuang (1982) predicts that the minimum stress intensity factor for steady-state crack growth is Kmjn = 0.83 M P a 6 which is roughly 1.7 times the critical intensity factor predicted by the Griffith theory. This stress intensity factor corresponds to a crack with a length of roughly 20 p m in this problem. It is of interest to study the crack evolution in this relatively low stress situation, because it is not something that can be inferred from Chuang’s steady-state results. Chuang (1982) found that the crack velocity is extremely sensitive to the stress intensity factor, and even when the stress intensity is only 20% above the minimum value, the crack is predicted to propagate the distance of 100 p m in just 30 minutes. The results presented here consider how a defect develops over a considerable period during the lifetime of a component. Figure 23a shows how the crack length increases over time. It can be seen that the crack grows at a stress intensity much less than that predicted by Chuang or the Griffith theory, albeit very slowly. The complete lifetime of the component in real time is approximately 16 years. From the initial development of the opening displacement in front of the crack in Figure 24a, one can see that matter is distributed over a large area for a small increase in crack length. Chuang’s analysis predicts that a steady-state crack traveling at its minimum velocity due to Kmin will plate material over a length of about 3 pm. In addition, such a crack will have a maximum opening displacement of roughly 5 nm at its tip and the crack width will decrease as the crack velocity increases. At first the crack travels very slowly and has a substantial width of about 100 nm. Therefore the increase in crack surface area is relatively small for the amount of material redistribution required

142

Alan C. I;: Cocks et a1

0.045

0 04 0 035

d(al 0.03 a0

0 025 0 02

0 015 0 01 0 005 0

1

1

5

10

t/to (b)

FIG.2 3 . The temporal evolution of(a) the crack length and (b) the crack tip opening displacement under a remotely applied load of 100 MPa.

Modeling Microstructure Evolution

143

to extend it. There is ample time for the material to redistribute over a significant fraction of the grain boundary and from Figure 24c it can be seen that this reduces the stress at the head of the crack very effectively. Figure 23b shows that, as the crack slowly grows, the opening displacement at the crack tip monotonically increases; i.e., the amount of material wedged into the boundary, and hence the total elastic stored energy, increases. By the time the crack reaches the critical length of 20 b m predicted for steady-state growth (at a speed that would cause the crack to completely propagate the length of the boundary in 20-30 days), substantially more material has been wedged in front of the crack (35 nm at the tip) than would be expected if it had time to adopt a steady-state distribution (5 nm at the tip). It can be seen that this effect due to the plating history of the boundary is quite significant. The large mass of material wedged in front of the crack presents a substantial hindrance to its forward motion. In fact, it can be seen from Figure 23b that material continues to accumulate in front of it for a while longer. Eventually, this reservoir of elastic stored energy loses it battle with the applied load and the crack finds a mechanism by which it can release this energy. Relatively suddenly, the crack becomes very thin and moves forward through the bulk of the material wedged in front of it, leaving i t behind as indicated by the sharp drop in the opening displacement in Figure 23b which is associated with the rapid increase in the crack growth rate seen in Figure 23a when the crack is about 30 p m in length. Figure 24d shows that a large peak in the stress develops at the crack tip and it can be seen from Figure 24c that a substantial amount of material is removed from the grain boundary to the crack tip in an effort to try to reduce this (this was also found by Chuang ( 1982) for high-velocity cracks). As the stress increases the opening displacement and stress become more and more localized in the vicinity of the crack tip until the component fails. Similar behavior is observed for greater applied loads and the length at which the crack growth rate rapidly increases is found to be approximately inversely proportional to this load. Here we have examined a relatively simple loading situation and cracked geometry. More complex loading situations and geometries of the cracked body are examined more fully by Gill ef ul. (1998). The analysis described in this subsection takes into account the influence of the elastic deformation on the grain boundary deposition process, but the effect of changes of elastic stored energy on the shape of the growing crack has not been taken into account. In the next subsection we describe how this can be modeled and examine a situation where the elastic stored energy has a significant effect on the way in which the microstructure evolves.

Alan C. E Cocks et al.

144

a=ao

- - - a=4ao

0.035 2:

-. -.-

a = llao ............. a=23ao

0.03- !: ,*

L

0.025 4. *: : ' 0.020.015

*'

-

0.01-

0.005

\*

.*.

\. *.'*... !, \

\

-

.\**-...

\;.* *.:\

">

\ \

\

01

- - _.._-._

'-

I

-.-. ........................................................... -.=-.=.= .-_-_-. -.-----

I

10

0

.:.>.

a ->.:

I

I

I

40

30

20

I

60

50

-a=40ao

0 02

- - - a=49ao

6(x) -

-._._ a=64ao

ao

............ a = 8 8 a o

0.01

t

.

"..(--"'................................................................................................ , -.- - -.-.-.-.-.- - -.- - - - - -.- - - - - - - - . __---*

\ *' I

-0.0'

c

*&

c

c \

-

-

-

_----

#

I

L'

-0.0,

I

1

I

2

I

3

I

4

I

5

I

I

I

I

6

7

8

9

10

FIG, 24. Illustration of (a)-(b) the opening displacement and (c)-(d) the slress field. ahead of a growing crack of finite length u under a remotely applied load of 100 MPa.

2.0

I

I

I

1

I

I

a=ao

1.8

- - - a = 4ao

- _ _ _a =_Ilao

1.6

............. a = 23ao

-

........................ ... ........... .... ..... .... -.-.-.-..... ............. -.-.....---._. -. .....-..... ...........

0 ( I:)1 .......... -

1..

-------

c---------

- - - - - ----_ -.

I' /

I'/

-

0.8 7

06

0

25 -;

I

I

I

1

I

I

10

20

30

40

50

60

I

I

I

I

I

$

20

4jl

I!

-! ! \ I

10

-a = 4 0 a 0

I

a = 49ao

-

a = 64a0 ............ a = 88ao

-

: I

\'I ~(x) * .

a

I

I

- - -

:I :. 15

I

70

**

.

- :. 1,

-

.\>..

.* \. ' Y. 5

5-

\*

**-*..-.

----.. '. - .- .-. - .y"-=............... - .................................. - _ -............ --

-

--A__

0-

I

I

I

1

I

I

I

I

I

146

Alan C. E Cocks et al. 2. Su f a c e Evolution of an Elastically Strained Epitaxial Thin Film

Thin films are widely used in high-technology applications, such as microelectronic and optoelectronic devices. They are generally manufactured by the deposition of a material from the liquid or vapor phase onto a supporting substrate. Under certain deposition conditions, epitaxial films are formed in which the atomic lattices of the film and substrate materials are commensurate. If no defects are formed, large elastic mismatch strains can result (roughly 4% for Si on Ge and up to nearly 7% for GaAs on TnAs) to ensure that the interface is coherent. The major contributions to the Gibbs free energy of the system are therefore the surface energy and the elastic stored energy. The elastic stored energy of the system can be reduced by the surface becoming wavy (Gao, 1991; Freund and Jonsdottir, 1993) but this results in an increase of surface energy. The driving force for the growth of these waves (which generally occurs by surface diffusion) and the ultimate breakup of the film (and hence failure of the device) is determined by the competition between these two contributions. We assume that the surface profile of the film evolves by surface diffusion, which is the only kinetic process. Therefore

The surface profile evolves in order to reduce its surface energy y 5 A , and the elastic strain energy due to the film-substrate mismatch strain U,, so that

G = y,A,

f

u,.

Assuming that the elastic strains are small, the elastic strain energy in the body can be determined through a transformation strain analysis in which the film is given a transformation strain E,; (the negative of the strain required to make the film coherent with the substrate). There are two contributions to A , : the surface distortion resulting from the elastic deformation of the film and substrate and changes of the surface profile arising from the diffusional redistribution of material. For the situations considered in this section, the contribution arising from the elastic deformation of the film is small. We are then able to describe the evolution of the profile of the film in terms of a single dimensionless parameter Eo, defined below. An evaluation of the effect of the elastic deformation of the body is given by Cocks and Gill ( I 998). For a film of volume V f on a substrate of volume V,, the elastic stored energy is given by

Modeling Microstructure Evolution

I47

where D;,kl is the elastic stiffness tensor and EL, is the total strain at a material point. Discretization of the film into a number of finite elements allows the elastic strain energy to be expressed as

(4.28) where all displacements and lengths are normalized with respect to the film thickness h , [ K,] is the dimensionless stiffness matrix for a given configuration, [ R ] is a normalized transformation strain matrix, [S]is the matrix of elastic displacements of the elemental nodes, and E , , ~ is the magnitude of the mismatch strain. Since we are only considering rearrangement due to surface diffusion, the last term of eq. (4.28) is constant throughout the evolution process. The nodal displacements are those that minimize the elastic strain energy:

I S 1 = ~ , , l K , l - ' l R l = &,,[GI. The rate of change of strain energy is then given by

where l u ] is the velocity matrix of the nodes on the surface of the film with normalized coordinates [x ], and driving force matrix [ F(,] determined from changes in the elastic strain energy is the assembly of

where F,, represents the contribution to the force matrix from the change in elastic stored energy associated with an increase of the ith degree of freedom. To evaluate this contribution, one must determine the differentials of the stiffness and transformation matrices. Analytical expressions for these differential matrices can be readily evaluated for isoparametric finite elements (Cocks and Gill, 1998)in terms of the mapping from the parent (reference) element to the mapped (physical) element. We normalize the variational functional in terms of the reference quantities

so that

I48

Alan C. F: Cocks et al.

where the dimensionless parameter Eo = Eh&:,/ y, represents the relative significance of the elastic mismatch strain energy and surface energy as driving forces. As in the analysis of surface diffusion in Section IV.B, we must modify this functional by introducing a series of Lagrange multipliers associated with the surface nodes to ensure that matter is conserved as the surface evolves. The nodal velocities of the film surface are given by the stationary solution of this modified functional:

and hence the surface profile of the film can be updated by a suitable time increment. Once this has been done, the finite-element mesh must also be reconstructed. Freund ( 1995) analytically investigated the evolution of small surface perturbations in an elastically strained thin film using a Fourier series representation for the surface profile, Results were obtained for a number of axisymmetric imperfections by the numerical integration of an infinite series of Bessel functions. The results presented here are for two-dimensional plane strain problems. Therefore, using the same notation as Freund. the perturbation from the uniform flat surface profile is chosen to be

A h ( x ) = -Aexp)Z: ; (

= -Aexp

)?I;(

__

,

where A and u define the shape of the imperfection and the asterisk represents the normalization x * = x / l for lengths and I*= r / t for time in which 1 = y , / U o = 2 h / E o is a reference length, Uo is the strain energy density in an unperturbed film, and t = l 3 / D 0 , yis , a reference time. Freund estimated that T x 10 s and 1 FZ 1 F m for silicon under a biaxial strain of 1% at 800 K. An initial imperfection of this type of amplitude 0.01h was introduced into a film of thickness h and width 2h and a suitable finite-element mesh was chosen. The results of the numerical simulation are shown in Figures 25 and 26 for two different imperfection sizes. Although these results in two dimensions are not directly comparable with those of Freund in three dimensions, their general form and behavior are very similar. The strain energy acts to increase the amplitude of the perturbation. However, when a* = 0.5 the imperfection is quite sharp, and the driving force to reduce the curvature at the center of the imperfection dominates initially, causing the amplitude to decrcase momentarily. This is accompanied by a spreading out of the imperfection along the length of the film. After this transient period the imperfection propagates toward the substrate. Using the method

Modeling Microstructure Evolution

I

t*=O ........... t*=0.05

I

-

-1.5

/

I

-.-.-

---

/

-2

4

-0.5y -2SL

I

1

0.5 1

I

I

-1

Ah A

149

/

I

I

I

1

I

I

1

I

I

I

t"z0.22 t*=0.43 I

I

-

I

I

,'..'

/r;

- 1.5

........4'

.' I

-t * =o

I

-2

............ t* =0.15 -.-.- t*=o.2

- - - t*=0.25 -3.5

0

I

0.5

I

1

I

1.5

I

2

I

2.5

I

3

I

3.5

I

4

I

4.5

X*

(b) F I G . 25. Surface profile evolution of a small perturbation for (a) (I* = 0.5 and (b) (I* = I .0

Alan C. E Cocks et al.

150

L bhlojl A -1.2

t ..

t

\

\ \

FIG.2 6 . Temporal evolution of the amplitude of a small surface perturbation.

presented here, the development of these perturbations into sizeable defects can be followed through to the breakup of the film or other geometrical configuration (Cocks and Gill, 1998). The procedures presented here can also be readily extended to consider the formation of continuous and discontinuous thin films by including the process of condensation from the vapor phase and accounting for the accompanying increase in volume of the film.

V. Rayleigh-Ritz Analyses In the previous section we described how numerical simulation procedures can be developed from the variational principle of Section 111. An accurate description of the geometry is often required to obtain a simulation that fully describes the way in which the microstructure evolves. We should bear in mind, however, that the geometric situations we consider in these simulations are necessarily idealized and when we compare the simulations with experiments we are often interested in global features of the evolution process. We recognized this in our modeling of grain growth in Section IV.F, where we focused on a limited number of features of the evolving microstructure and used the variational principle to develop constitutive laws in terms of a small number of state variables and global macroscopic parameters. Also, diffusion coefficients are often not known with any degree of

Modeling Microstructure Evolution

(a)

151

(b)

FIG.27. The cosintering of two spheres showing the detailed shape of the particles determined using the techniques described in Section 1V.B. (a) Early in the sintering process, showing the formation of a neck. (h) Late i n the process. when diffusional rearrangement of material over the entire surtace of the particles has occurred.

accuracy; surface diffusion coefficients can only be taken as order-of-magnitude estimates. In any modeling exercise we should determine the overall objectives of the modeling and determine the simplest description that retains the major features of the evolving microstructure. Simulations can then be undertaken in an efficient manner and larger-scale simulations can readily be undertaken which allows information to be passed up the full hierarchy of length scales. Here we demonstrate how relatively crude descriptions of the evolving microstructure can be employed to provide insight into the dominant physical processes. In Section 1V.B we examined the cosintering of a row of identical spherical particles. The evolving shape of two contacting spheres computed using the methods described in Section I1 is shown in Figure 27. Early in the sintering process a neck region develops (Figure 27a). Sintering occurs by material flowing away from the grain boundary formed between two adjoining particles and plating out on the surface of the particles in the vicinity of the neck, with the detailed shape of the neck region determined by the surface and boundary diffusion processes. Away

Alan C. E Cocks et al.

I52

Frc;. 2 8 . An idealized representation of Figure 27. with the geonietry represented by truncated spheres of radius R separated by cylindrical disks of thickness t and radius .r.

from the neck region, the particles retain their initial equilibrium spherical shape. As sintering progresses the neck region grows and material rearrangement occurs over the entire particle surface (Figure 27b). Eventually, an equilibrium configuration develops in which the surface has a constant curvature. The early stage of sintering is similar to that experienced by a row of contacting spheres. Parhami er al. (1998) modeled this situation by adopting a simple two-degree-of-freedom geometric description (Figure 28). At a given instant the geometry is represented by a row of truncated spheres of radius R separated by disks of radius x and thickness d , with a grain boundary located at the center of the disk. The condition that volume must be conserved provides a relationship between x , d , R , and the initial radius Ro. The disk represents the development of the neck region and changes in the radius and thickness of this region indicate the extent of the local diffusion process. Changes in the value of R reflect the relative importance of the global diffusion process. Parhami et al. (1998) computed the evolution of the simple profile using the variational principle of (4.2) and compared the resulting evolution of the neck radius, x, and center-to-center spacing with more detailed in the range 0.1 to 10 for both free and forced sintering. These calculations for general features of the evolving microstructure are in very good agreement with the more detailed calculations over the full evolution process. Here we consider in more detail the early stages of the sintering process under the action of an axial force F . During the early stages of sintering, material only diffuses local to the neck region. As a result, as observed by Parhami et al. (1998), R does not change and remains equal to its initial value Ro. Thus we can treat the problem as a singledegree-of-freedom problem. Volume conservation gives, for x and d much less

es

Modeling Microstructure Evolution

153

than Ro,

We express the rate of change of the microstructure in terms of the height h of the unit cell shown in Figure 28, where 7

7

h=2[R,j-x-]

112

+d=2Ro-d

(5.2)

and = -2.The flux along the circular grain boundary can be obtained from the continuity condition of (3.3) and the flux along the surface of the cylinder from the continuity condition

Substituting the resulting flux fields into the variational principle of (3.6), noting that

1 . A/, = - - A , = 2 2

~

~

i

for small x and d and optimizing with respect to 6, gives

or

In general, 3,is greater than 1 and for the conditions considered here x' << Ro. The term involving 3,can therefore be ignored, and integration of (5.3) from x' = 0 to its current value gives

(5.4) For F = 0 and y/, = 0, eq. (5.4)reduces to Coble's (1958) free sintering result. We can therefore regard (5.4) as a generalization of Coble's analysis. Parhami et al. (1998) demonstrate that this expression is in good agreement with their fuller simulations, in which they also allowed R to vary, and more complex numerical models up to values of x / Ro of the order of 0.1. For larger values of x / R o , diffusion over the entire particle surface becomes significant and the rate of neck growth decreases as the geometry approaches the equilibrium configuration.

154

Alan C. E Cocks et al.

The real power of this approach is not that we are able to model this problem with a reasonable degree of accuracy using only a small number of state variables, but it allows us to model problems over a much wider range of length scales, i.e., to analyze large arrays of particles. This does not mean that there is no role for detailed computations of the type described in Section IV. It is only by performing simulations of that type that we are able to identify the major features of the way in which the microstructure evolves. Thus simple models of the type described here are guided by the physical insight developed from the more detailed simulations. Parhami (1 996) has used the simple model described here as a major building block in the development of a discrete element model for the sintering of random arrays of particles. In the following section we use the result developed here to obtain an analytical constitutive model of the sintering of ceramic compacts.

VI. Structure of Constitutive Laws for the Deformation of Engineering Materials In Section IV we used the variational principle of Section 111 to develop numerical procedures for the detailed modeling of microstructure evolution in engineering materials. In many practical situations we wish to develop constitutive laws for the behavior of a material in terms of a limited number of state variables. In this section we identify a general structure for constitutive laws where inelastic deformation is due to the class of diffusional mechanisms considered in Section 11. In the process we simplify the structure of constitutive laws proposed by Cocks (1994). Consider a body of unit volume which is subjected to macroscopic stresses Xi,, and experiences inelastic straining at a rate E i j . We limit our discussion to situations where the free energy is dominated by the potential energy of the applied load and the energy associated with any interfaces in the body. For this situation the rate of change of Gibbs free energy is given by

For situations where the interfaces consist of grain boundaries and surfaces and the specific energies are independent of state, we can simplify this relationship:

Modeling Microstructure Evolution

155

Substituting for G in eq. (3.5) and differentiating the functional to obtain the exact solution gives Cij dEij

- YbdAb - y.5 d A s

=d q .

(6.2)

In general, the interface areas are a function of the inelastic strain. Then

where Xi’/is a capillarity, or sintering, stress, which arises from the change in interface area that accompanies an increment of strain, i.e.,

+

where A = ybAb ysA,. Similarly, substituting (6.1) into (3.7) and considering chemical potential and stress fields that satisfy the surface and boundary energy requirements, we obtain

E,, dC,, =d @ , i.e.,

(6.5)

a@ 8..= IJ axij

’

The rate of energy dissipation per unit volume is given by

D

= -G = @ + P ‘

and the macroscopic work rate per unit volume is lQ = X,j E;, = 0

+ \v + A .

These results simplify the expressions originally proposed by Cocks (1994) by relating all contributions to the capillarity stress to changes in internal interface areas and combining all these contributions into the scalar potential A. Following Cocks (1994), we can obtain bounds on the macroscopic potentials. We do not repeat the complete details here and further limit our attention to the kinematic bounds derived from the variational principle of (3.6). Bounds on 0 and Q A can be obtained by identifying suitable compatible flux and velocity fields, which results in a macroscopic strain rate E c . For a body subjected to a prescribed stress C i j , the strain rate potential is bounded from below by

+

@

>_ & j E Y .‘ J - Q c - A‘,

(6.7)

Alan C. F: Cocks et al.

I56

where the superscript c indicates quantities determined from the assumed compatible field. If instead the strain rate k,, i5 prescribed and the assuined field is scaled such that hi, = E , , , then

W

+ A 5 W' + A ' .

(6.8)

These bounds and their extensions allow the nature of the interactions between different mechanism$ to be determined and suitable simplified structures for constitutive laws to be identified. Full details are given by Cocks (1994). In this section we demonstrate the utility of the above bounds by developing an anisotropic constitutive law for stage 1 sintering controlled by grain boundary diffusion employing the bound of (6.8).

A. STAGE1 COMPACTION OF CERAMIC COMPONENTS Consider the random array of spherical particles of initial radius Ro shown in Figure 29a, which has been strained to a state E,, and is subjected to a strain rate E , , . We limit our attention to the early stages of compaction, when the strains and contact patches formed between the particles are small. We further restrict our consideration to conditions where all the principal components of strain, but not necessarily strain rate, are compressive. It was demonstrated in Section V that, if the necks formed between the contacting particles are small, diffusional processes only occur local to the neck region. We can therefore assume that there are limited interactions between the necks. Following the analysis of Section V, we idealize

(a)

(b)

FIG. 29. (a) A random array of spherical particles subjected to an applied strain rate conlacring particles from [he array of (a) with outward normal n.

b;,.(h) Two

Modeling Microstructure Evolution

157

the geometry at a given instant as an array of truncated spheres of radius Ro separated by disks of material of radius x and thickness d . We assume, as in the analysis of Fleck et al. (1997),that we can relate the center-to-center spacing of two particles whose contact patch has an outward normal n; (Figure 29b) to the strain: Nn,)=2Rdl

+ E,,n,n,).

(6.9)

During the early stages of compaction, the average number of contacts per particle, Z, will remain essentially constant. Geometric considerations allow us to relate the contact radius x ( n , ) to h ( n , ) (see Section V) and hence E l , : x ( n , ) = 2Ro(-E,JnlnJ)'/'.

(6.10)

We now assume a compatible velocity field such that h(n,')= 2R&,,nln,

(6.1 1)

and, from the analysis of Section V, the contribution of the boundary diffusion process on the contact patch to the potential W is (6.12)

Consider an element of the compact which contains a single particle and has a volume :rr R i / D , where D is the relative density of the compact. Since the array of particles is random, the number of contacts per unit area is

z=- Z

4rr R,'

and

(6.13)

W becomes (6.14)

where the factor of arises because a given boundary is shared between two particles and S is the total surface area of a particle. Substituting for z and @'(ni) using (6.12) and (6.13) gives (6.15)

To determine the capillarity stress contribution, we need to relate the rate of change of the area of the contacts to the applied strain rate. Employing the above

158

Alan C. F Cocks et al. COMPONENTS O F VISCOSITY

TABLE2 M A T R I XDEFINIIIII N ( 6 .18): cj,,i\=

670K +c,,ji[

compatible field and making use of the analysis of Section V, we find 1 . 7 ‘ A b ( n , )= - - A s h f ) = -4rrR6El,n,n, 2

(6.16)

and

where EL,is the volumetric strain rate. Choosing a coordinate system whose axes lie in the principal strain directions, combining (6.15) and (6.17), and evaluating the integrals of (6.15), we find

where

The viscosity matrix CfJx/can be expressed in terms of the three principal strain components and four constant parameters to give six independent components of the viscosity matrix. All the nonzero components are given in Table 2. It is interesting to note that the sintering potential defined in (6.18) is only a function of the volumetric strain (i.e., relative density of the component). This does not, however, mean that the material will densify uniformly after experiencing a prior loading history which results in an anisotropic structure. The vis-

Modeling Microstructure Evolution

159

cosity in the direction of the maximum compressive stress is greater than in the transverse directions. As a result, the material will strain at a slower rate in this direction and the material will deform toward a more isotropic state. If the material is compressed isostatically, the three strain components are equal and the material response is isotropic. The macroscopic response can then be expressed in terms of a single state variable, the relative density D . The general form of the model then reduces to that developed by McMeeking and Kuhn (1992). The approach presented here can be readily extended to other mechanisms of material redistribution. For isolated contacts, once the solution for the deformation between two contacting particles is known, the strategy described here can be employed to determine the macroscopic response.

VII. Concluding Remarks In this paper we have described a general variational principle for the analysis of microstructure evolution in engineering materials. We have demonstrated how self-consistent numerical schemes can be developed from it which naturally take into account the competition and synergy between the different mechanisms. The form of the variational principle allows the relative importance of different mechanisms to be readily identified and strategies to be developed in which the full range of possible situations can be systematically evaluated. In Section IV we limited our attention to situations in which material rearrangement occurs by direct diffusional processes or through the migration of interfaces. In the process we have considered a wide range of thermodynamic driving forces. The full range of kinetic processes and the origin of the thermodynamic driving force are summarized in Table 1. This is not an exclusive list. The techniques described here can be extended to a wide range of mechanisms. For example, power law viscous flow of the matrix and multicomponent diffusion problems can be considered within the current framework as well as other contributions to the driving force, such as chemical contributions related to composition and phase changes. Similarly, the range of physical problems we have examined is not intended to be exhaustive. These have been chosen to illustrate the techniques and to demonstrate how macroscopic constitutive laws can be developed in a systematic manner from the results of the simulations. We have also shown how simplified methods of analysis of the Rayleigh-Ritz type can be developed. These provide a good indication of the general way in which the microstructure evolves. Care, however, needs to be taken in the use of these methods. Enough degrees of freedom need to be chosen to allow for the different major ways in which the microstructure can evolve. Use of this method

160

Alan C. E Cocks et al.

therefore generally requires a deeper physical understanding of the way in which a given material system is likely to behave than use of the full computational methods. It has also been shown how bounding theorems can be developed from the variational principle of Section 111. A major use of these bounds has been in the development of constitutive laws for the inelastic behavior of engineering materials, as described in Section VI. A fuller evaluation of the use of these bounding theorems is provided by Cocks (1994).

References Amerasekera. E. A,, and Nam, F. N. (1997). Failure Mechanisms in Semiconductor Devices. Wiley, Chichester. Arzt. E., Ashby. M. F., and Verrall, R. A. (1983). Interface controlled diffusional creep. Acra Mefall. 31. 1977. Ashby, M. F. (1969). On interface reaction control of Narharro-Herring creep and sintering. Scripra Meftrll. 3, 837. Ashby, M. F., and Verrall, R. A. (1973). Diffusion accommodated flow and superplasticity. Acra Metall. 21, 149. Atkinson. H. V. (1988). Theories of normal grain growth in pure single phase systems. Acra Mefall. 36, 469. Burke, J. E. (1949). Some factors affecting the rate of grain growth in metals. Trans AIME 180.73. Burton, B. (1972). Interface reaction controlled diffusional creep: a consideration of grain-boundary dislocation climb sources. Marer: Sci. Engrg. 10, 9. Chen, I. W., and Xue, L. A. (1990). Development of superplastic structural ceramics. J. Amer: Ceram. Soc. 73. 2585. Chuang, T.-J. (1982). A diffusive crack growth model for creep fracture. J. Amer: Cerani. Soc. 65. 93-103. Chuang, T.-J., Kagawa, K. I., Rice, J. R., and Bank-Sills, L. (1979). Non-equilibrium models for diffusive cavitation of grain interfaces. Acra Mefull. 27, 265-284. Chuang, T.-J., and Rice, J. R. (1981). Energy variations in diffusive cavity growth. J. Amer: Ceram. Soc. 64, 46-53. Coble, R. L. (1958). Initial sintering of alumina and hematite. J. Amer: Ceram. Soc. 41, 55. Cocks, A. C. F. (1990). A finite element description of grain-boundary diffusion controlled processes in ceramic materials. In Applied Solid Mechu,iics--3 (I. M . Allison and C. Ruiz, eds.). Elsevier, North-Holland, Amsterdam. Cocks, A. C. F. (1992). Interface reaction controlled creep. Mech. Mater: 13, 165. Cocks. A. C. F. (1994). The stmcture of constitutive laws for the sintering of fine-grained materials. Actcr Metall. Mate,: 42, 2191. Cocks, A. C. F. (1996). Variational principles, numerical schemes and bounding theorems for deformation by Nabarro-Herring creep. .I Mech. . Phys. Solids 44, 1429. Cocks, A. C. F., and Ashby. M. F. (1982). On creep fracture by void growth. Prog. Mater: Sci. 27, 189. Cocks. A. C. F., and Gill. S. P. A. (1996). A variational approach to two dimensional grain growth. 1. Theory. Acta Murer: 44. 47654775. Cocks, A. C. F., and Gill, S. P. A. (1998). A numerical model for stress-driven diffusion in elastically strained epitaxial thin films. Leicester University Engineering Department Report 98-3.

Modeling Microstructure Evolution

161

Cocks. A. C. F.. and Pan. J. i 1994). Thc influelice o f an interface reaction on the creep response daniaged materials. M t 4 Mrrwr: 18. 269. Cocks, A C . F.. and Searlc. A. A. 11990). Cavity growth in ceramic materials tinder multiaxial stress states. A C I ~MJ r ~ l l38. . 2493. Cocks. A. C. F.. and Searle, A. A. 119911. Void growth by grain-boundary diffusion i n tine grained materials. Meelr. M ~ J I 12, ~ J 279. : Doheny. R. D. 11975). Discusion of Mechanism of steady-stair grain growth i n Aluminium. Metrrl-

/JlJ:qfCrf/ Trtrtrv. 6A. 588-590. Du. Z.-2.. and Cocks. A. C. 1;. (1992). Constitutive models for the sintering ofcet-amic components. 1. MLiterial models. Acm Mmrll. 40. 1969. Du. 2.-2..McMceking. R. M.. and Cocks. A. C. F. i1998). A numerical model of grain growth. Philos. Mtrg. To appear. Fleck. N. A,. Storikers. B.. and McMeeking. R. M. ( 1997). The viscoplastic compaction of powders. 111Mrchtrr~icsof Grtrrirrkrr t r ~ i r Po,nir.\ l Mrrtt~ricrls(N. A . Fleck and A. C. E Cocks, eds.). Kluwer Academic, Dordrecht. Floro. J. A,. and Thompson. C. V. i 1993).Numerical analysis o f interface energy-driven coarsening in thin films and its connection to grain growth. Arm Metnll. Mtrtrr: 41. 1137-1 147. Freund. L. B. i 1995). Evolution of waviness on the suilace of a strained elastic solid due to stressdriven surface diffusion. I r m w i r r / . J. Solid.\ Srrrrcnrres 32, 9 1 1-923. Freund, L. B., and Jonsdottir. F. i1993).Instability o f a hiaxially stressed thin tilm on a substrate due to materials diffusion over its free surface. J. Mecli. Phw. Solir/.s 41, 1245-1264. Gao. H. ( 1991). A boundary perttirhation analysis for elastic inclusion\ and interfaces. Irrternot. J. Solirl.~Strrrcrrrres 28. 703-725. Gluier. J. A.. and Weaire. D. (1992). The kinetics of cellular patterns. J. P l y Co!it/ert.st,d Mrrftrr 4, 1867. Gill. S . P. A,. and Cocks. A. C. E i1996). A variational approach to two dimensional grain growth. 11. Numerical results. Ac/tr Mrrwc 4 4 . 1 7 7 7 1 7 8 9 . Gill. S. P. A,. and Cocks. A. C. F. (1997a). An investigation of mean-field theories for normal grain growth. Plrilos. Mtrg. A 75. 30 1-3 13. Gill, S. P. A,, and Cocks. A. C. F. (1997b). A short note on a variational method for two-dimensional normal grain growth. Scrip~trMLJWJ: 35. 9. Gill. S. P. A,. and Cocks. A. C . F. ( 1998). A variational approach to modelling abnormal grain growth i n thin tilnis. Leicester University Engineering Report 98- I. P. A.. Cocks. A. C. F.. and Gao, H. (1998). A variational model of the diffusive growth of an k of finite length. Leicester University Engineering Department Report 9X-2 Heywang. W. (197 I). Semiconducting barium titanate. J. Mtrtrt: S r i . 6, 12 14. Hillert. M. 11965). A mean field theory foi- nornial grain growth. Amr Mrtrill. 13. 227. Honeyconihe. R. W. K., and Bhadeshia. H. K. D. H. 1 1995). Sfrrls~Mi~.ro.strrrc/rrrr mid Proprrrirs. Edward Arnold. London. Howes. M. J.. and Morgan. D. V. ( 19851. Gtrllirtt~rAtwrrielr: Mrr/t,riu/s. Drvicrc rirttl Cinwits. Wiley. Chichester. Kucherenko. S.. and Pan. J. 1 1998).A coupled linite element and finitc diffcrcncc formulation with an implicit time integration scheme for analysih of inicrwtructure evolution. LJniversity of Surrey Report. Luo, K. K . (1978). Analysis of branched cracks. J. Appl. Mrch. 45. 797-802. Matan. N.. Winand. H. M. A,, Bogdanufl. P. D.. and Reed. R. C . (1998). Coupled thermodynaniickinctic modelling of diffusion reactions i n superalloys, Actrr Mnrrr. To appear. McMeeking. R. M.. and Kuhn. L. T. i1997).A diffusional creep law for powder compacts. Acrrr M e ~ l l . Mrrtc.~:40. 96 I.

162

Alan C. E Cocks et al.

Mullins, W. W. (1958). The effect of thcrmal grooving on grain-boundary migration. Acfa M ~ t u l l6, . 414. Pan, J., and Cocks, A. C. F. (1993a). The effect of grain-size on the stress field ahead of a crack in a material which deforms by Coble creep. Internut. J. Fracture 60, 121. Pan. J., and Cocks, A. C. F. (1993b).Computer simulation of superplastic deformation. Compuf.Mare): sci. 1,95. Pan, J . , and Cocks, A. C. F. (1995). A numerical technique for the analysis of coupled surface and grain-boundary diffusion. Act0 Metall. Muter: 43, 1395. Pan, J . . and Cocks, A. C. F. (1998). Modeling crack growth in engineering ceramics. To appear. Pan, J . , Cocks, A. C. F., and Kucherenko, S. (1997). Finite element formulation of coupled grainboundary and surface diffusion with grain-boundary migration. Proc. Roy. Soc. London Se): A 453,2161. Parhami, Z. (1996). Ph.D. thesis. University of California at Santa Barbara. Parhami, F., McMeeking, R. M., Cocks, A. C. F.. and SUO,Z. (1998). A model for the sintering and coarsening of rows of spherical particles. Submitted for publication. Perduijn, D. J., and Peloschek, H. P. (1968). Mn-Zn ferrites with very high permeabilities. Proc. British Ceram. Soc. 10, 263. Pharr, G . M., and Nix, W. D. (1979). A numerical study of cavity growth controlled by surface diffusion. Actu Mefall. 27, 1615-1631. Raj. R., and Ashby, M. F. (1975). Intergranular fracture at elevated temperature. A c f a Mefall. 23,653. Raj. R., Shih, H. M., and Johnson, H. H. (1977). Correction to: Intergranular fracture at elevated temperature. Scripta Merull. 11, 839. Rhines, F. N., and Craig, K. R. (1974). Mechanism of steady-state grain growth in Aluminium. Merullurgical Trans. 5.413425. Shewmon, P. G . (1963). Difusion in Solids. McGraw-Hill, New York. Spriggs, R. M., and Dutta, S. K. (1974). Pressure sintering-recent advances in mechanisms and technology. Sci. Sintering 6 , 1. Sun, B., SUO,Z . , and Cocks, A. C. E (1996). A global analysis of structural evolution in a row of grains. J. Mech. Phys. Solids. Suo, Z. (1996). Motions of microscopic surfaces in materials. Adv. Appl. Mech. 33. Thompson, C. V. (1985). Secondary grain growth in thin films of semiconductors: theoretical aspects. J. Appl. Phys. 58,763-772. Thompson, C. V. (1990). Grain growth in thin films. Ann. Rev. Mates Sci. 20, 245-268. Thouless. M. D., Hsueh, C. H., and Evans, A. G . (1983). A damage model of creep crack growth in polycrystals. Aclu Memll. 31, 1675-1 687. Van Siclen, C. Dew. (1996). Random nucleation and growth kinetics. Phys. Rev. B 54, 11845-1 1848.