A computer simulation method for particle sintering

Acta metall. Vol. 36, No. 8, pp. 1977-1987, 1988 Printed in Great Britain. All rights reserved

0001-6160/88 $3.00 +0.00 Copyright ~ 1988 Pergamon Press pie

A COMPUTER SIMULATION METHOD FOR PARTICLE SINTERING H. J. LEU, T. HARE and R. O. S C A T T E R G O O D Department of Materials Science and Engineering, North Carolina State University, Raleigh, NC 27695-7907, U.S.A. (Recewed29July 1987; m rev~edform 11 January 1988) Abstract--A new computer simulation method is developed for treating topological effects during early stage sintering in large-scale particle arrays. The approach is based on a stress-free transformation concept using two-body sintering models where elastic distortions are introduced in order to retain the particle network connectedness. The internal stresses associated with these distortions can drive creep mechanisms. An initial implementation of the simulation method has been made for two-dimensional particle arrays, but it can also be applied in three dimensions. Empirical parameters were used to demonstrate the viability of the method. The results showed that the elastic stiffness parameters of the particle network can affect topological rearrangement as well as the distribution of internal stress. A sintering-creep interaction was found such that, when creep effects are present, the densification rate can be significantly reduced. R6sum6----Nous d6veloppons une nouvelle m6thode de simulation sur ordinateur, pour traiter les effets topologiques pendant les premiers stades du frittage darts des arrangements de particules ~i grande 6chelle. Cette approche est bas6e sur un concept de transformation sans contrainte, qui utilise des moddles de frittage ~i deux corps oft l'on introduit des distorsions 61astiques pour maintenir la cohesion du r~seau de particules. Les contraintes internes associ6es ~i ces distorsions peuvent conduire fi des m6canismes de fluage. Une premi6re mise en oeuvre de la m6thode de simulation a 6t6 faite pour des arrangements bidimensionnels de particules, mais on peut l'appliquer 6galement au cas de trois dimensions. Nous avons utilis6 des paramdtres empiriques pour d6montrer la validit6 de cette m6thode. Les r6sultats montrent que les constantes ~lastiques du r6seau de particules peuvent avoir un effet sur les r6arrangements topologiques tout autant que sur la r6partition de la contrainte interne. Nous avons observ6 uric interaction fluage-frittage capable, en pr6sence d'effets de fluage, de r6duire nettement la vitesse de densification.

Zusammenfassung--Es wird ein neues Computerverfahren zur Simulation topologischer Effekte entwickelt, welche im Friihstadium des Sinterns in groflr~iumigen Teilchenanordnungen auftreten. Das Verfahren benutzt das Konzept einer spannungsfreien Transformation auf der Grundlage von Zweik6rper-Sintermodellen, bei der elastische Verzerrungen eingefiihrt werden, um die Verbindungen im Teilchennetzwerk aufrecht zu erhalten. Die mit diesen Verzerrungen zusammenh/ingenden inneren Spannungen k6nnen treibende Kr/ifte fiir Kriechmechanismen darstellen. Eine erste Anwendung der Simulationsmethode betraf zweidimensionale Teilchenanordnungen; sie kann aber auch auf dreidimensionale angewendet werden. Anhand empirischer Parameter wurde die Richtigkeit des Verfahrens nachgewiesen. Aus den Ergebnissen folgt, dab die Parameter der elastischen Steifigkeit des Teilchennetzwerkes sowohl die topologische Umlagerung als auch die Verteilung der inneren Spannungen beeinflussen k6nnen. Es wurde eine Wechselwirkung zwischen Sintern und Kriechen gefunden, die bei Vorliegen von Kriechen zu einer betr/ichtlich verringerten Verdichtungsrate ffihren kann.

INTRODUCTION Theoretical analyses for early stage sintering (Stage I) have, for the most part, focussed on two-body sintering models without considering the important effects due to topological constraints in large-scale, random particle arrays. Topological constraints can cause particle rearrangements that influence the development of the void distribution in later stages o f sintering. Topological constraints can also lead to the development of internal stresses within the array. These stresses provide the driving force for creep mechanisms, which interact with the basic sintering mechanisms. Because of the inherent complexity and statistical nature o f the various phenomena involved, it is unlikely that any simplified analytical approach

can be developed to successfully deal with topological effects. The primary focus o f the work reported here was the development and initial implementation of a new, physically-based computer simulation method to treat topological effects during early stage sintering. Figure ! shows the idealized geometry used for two-body sintering models. Many such models have been reported in the literature, incorporating a variety of mass transport mechanisms to produce densification [1]. While such two-body models provide the basic kinetics for sintering, the sintering response of a large-scale particle array must differ from that for the two-body model. As will be discussed further in the next section, a given pair of particles in a large array of sintering particles will not

1977

1978

LEU et al.: COMPUTER SIMULATION OF PARTICLE SINTERING

FC F (a)

(b) Fig. 1. Two-body sintering model. The initial link length 1 is reduced by Al during a sintering time step At. Material forming the neck region is not shown. The topological constraint forces F do not affect the intrinsic sintering response; as explained in the text. These forces can be assumed to produce axial creep response superposed with the sintering response.

be free to undergo the same unconstrained shrinkage motion that occurs for an isolated particle pair. Topological constraints due to surrounding particles can lead to particle rearrangement and the buildup of internal forces or moments on the particles. Topological constraint effects in the early stage of sintering have been reported by numerous investigators [2-7]. Included among these are several important studies by Exner and co-workers [2-5]. These investigations have clearly shown that particle rearrangements occur during Stage I sintering, and that certain regions of a random, large-scale particle array will densify at the expense of other regions. The build up of internal stress within a sintering array of particles has been inferred from the fact that necks formed between particle pairs can be observed to distort and even fracture as sintering progresses [8]. Various computer simulation models for sintering have been reported in the literature [9-15]. The only previous computer simulation study that deals specifically with topological effects in large-scale randomly packed particle arrays was due to Ross et al. [9]. These authors clearly demonstrated that computer simulation could reproduce particle rearrangement qualitatively similar to that observed experimentally [2-5]. They also found that the initial particle packing arrangements can influence the sintering behavior. However, their simulation method was empirical in nature and effects such as topologically induced stress could not be treated. The fact that internal stress can be generated by topological constraint during sintering has been recognized in several recent studies [16-19]. DeJonghe et al. [! 6] and Hsueh [19] demonstrated experimentally that the effects can be very significant for the sintering of composite powder mixtures. Raj and Bordia [20]

and Hseuh et al. [17-19] attempted to model the stress effects analytically. Recently, Lange [15] developed an interesting analytical model that specifically recognizes the statistical nature of the internal stress due to topological constraint in a particle array. All existing analytical treatments are, however, highly simplified in their approach. A new computer simulation method for treating topological constraint effects in large-scale particle arrays during Stage I sintering will be discussed in the next section. Two-dimensional particle arrays will be treated because the concepts are most transparent for this case, however, the method is generally applicable to three-dimensional arrays. An initial implementation of a computer simulation code has been made to demonstrate the viability and generality of the method. SIMULATION METHOD Part&le arrays--network

representation

Random, two-dimensional particle arrays made up of spheres of uniform diameter shown in Fig. 2(a) were used for the simulation. The geometry and topological aspects of the array can be most easily

(a)

12 3

(b) Fig. 2. Representation of a randomly packed particle array (a). Network corresponding to the array (b). The link lengths and vertex angles are denoted by {lj) and {0j} as shown. The labeling scheme is arbitrary and for clarity not all links and angles are labeled.

LEU et al.:

COMPUTER SIMULATION OF PARTICLE SINTERING

1079

densification kinetics into the simulation code. The particular model used gives the two-body sintering law AI

r = A t t = 12 exp(l)

Fig. 3. Representation of bridging links in a network. Removing the bridging link separates the network into two distinct connected components. Two bridging links are shown (unfilled vertices) and the smaller connected component for each is identified by the dashed particles. A bridging link can connect to a single particle as shown. A sequence of bridging links terminating on a single particle forms a chain. represented using the network correspondent shown in Fig. 2(b). Particle centers are identified by vertices and a link joins a pair of vertices if the two particles are in contact. The network is thus made up of a unique, interconnected set of closed polygons. The detailed geometry of the network is given by the set of link lengths {lj} and angles {Oj} shown in Fig. 2(b). A map or table of the various interconnections between links and vertices specifies the topological properties of the network. Two different conditions on the outer array boundaries are possible. Free boundaries can be used, in which case free surface effects will be included in the simulation results. Alternatively, with suitable allowance for shrinkage effects, periodic boundary conditions can be used. The latter should better represent the behavior of an infinite particle array for simulations made with a limited number of particles. For computational convenience, two restrictions are imposed on the topology of the particle arrays. First, the network shown in Fig. 2 must be fully connected in the sense that there is at least one sequence of links connecting any two vertices. Second, the network should not contain bridging links (or chains) between connected components of the kind illustrated in Fig. 3. Loosely packed random arrays of particles contain bridging links while more densely packed arrays are relatively free of them. Because of these restrictions, the initial implementation of the method is restricted to more densely packed arrays. Such restrictions are not, in general, necessary. Sintering models, distortion

topological constraint and elastic

An empirical model [9] was chosen to introduce

(1)

where AI is the decrease in particle-pair spacing I for a given time step At (Fig. 1). Equation (1) is strictly an empirical form without consideration of physical parameters, and the choice of time or length scales is therefore arbitrary. An important feature of the simulation method rests in the fact that the sintering law can be arbitrarily chosen. The sintering simulation is based on the concept that a two-body sintering law can be applied to each contacting particle pair in the particle array. For a given time step At, each link/j in the particle network in Fig. 2(b) is subject to the sintering transformation l j = lj - Alj where As is determined using equation (1). This same concept was used by Ross et al. [9]. A central problem in developing the sintering simulation involves topological constraint. The physical origin of the constraint problem can be understood by considering the small cluster of four interpenetrating particles shown in Fig. 4. The network correspondent for the cluster in (i) is initially connected. To execute a sintering time step At, each link length in (i) is shortened using the two-body sintering law to generate the new link lengths in (ii). In general, the shortened links cannot remain topologically connected under the transformation so that initially closed and connected polygons in (i) will possess closure errors shown in (iii). Such closure errors produce the essential topological constraint. Closure errors shown in Fig. 4 will be accommodated in real particle arrays by elastic distortion

i

ii

iii

Fig. 4. A four particle network illustrating closure errors due to topological constraints. Polygon links in i are transformed during a sintering time step to those shown in ii. The orientation of the links is held fixed. If the links in ii are reconnected with fixed orientation, the closure errors in iii result. The initial interconnection of links in i can be retained in iii only if the links and vertex angles are elastically distorted.

LEU et al.: COMPUTER SIMULATION OF PARTICLE SINTERING

1980 (

~

F

Au

F

(a)

M (b) Fig. 5. Representation of axial (a) and bending (b) distortions for network links and vertex angles, respectively. The elastic stiffness coefficients k and//relate the distortions Au and Aco, respectively, to the internal forces F and moments M as explained in the text.

(stress). The shortened links in (iii) can be constrained so as to retain the original closed and connected polygon topology if elastic distortions occur during the sintering time step. That is, by elastically stretching, compressing or rotating the shortened links in (iii), the properly connected polygon topology in (i) can be restored, albeit with some change in polygon shape along with a concomitant amount of sinteringinduced shrinkage. A two-body sintering model is thus viewed as a stress-free transformation strain applied to a particle network. Connectedness of the network imposes a compatibility condition that produces elastic distortion. The concept is similar to transformation strain techniques used in continuum mechanics [21]. The mathematical formulation of the simulation method requires that linear elastic stiffness parameters be included into the physical model. As shown in Fig. 5, the axial stiffness coefficient k is defined such that if a link undergoes axial distortion Au, then the force F generated in the link is kAu and the stored elastic energy is l/2kAu 2. Similarly, if an adjacent pair of links joined to a common vertex undergo bending (rotational) distortion Ato, then the bending stiffness coefficient fl is defined such that the moment M generated in the common vertex is flAtn and the stored elastic energy is 1/23Ato 2. The inclusion of a proper treatment of vertex angle change is an important feature of the current approach. This has not been addressed in previous simulation methods [9]. The determination of actual numerical values for the stiffness coefficients requires detailed analysis of the elastic response of interpenetrating particles with given geometry and elastic properties, for example, using finite element techniques. As particle interpenetration increases during sintering, one can expect

that the stiffnesses will increase. Such detail was not attempted for the work reported here, and constant values of k and fl will be used for illustrative purposes. The elastic distortions generated in the particle network as a result of topological constraint are evaluated by making the fundamental assumption that an array of sintering particles remains in quasistatic equilibrium. This means that the inertial forces due to particle motion are negligible (zero inertial mass). The equilibrium distribution of the elastic distortions can then be determined by minimization of the stored elastic energy U. Letting {Auj} and {Atoj} represent the axial and bending distortions needed to maintain topological connectedness of the particle network at any time during the sintering history, U is given by 1

v

1

=-E Em + ,J

2j

(2)

where the first and second summations extend over all network links ~ and vertex angles 0r respectively. The reference state from which the elastic distortions {Auj} and {Ao~j} are measured is determined by the stress-free sintering transformation. The reference length for a given link in any current state of the particle network is the length obtained by subtracting from the initial link length all increments Al given by equation (1) for all sintering time steps At that have been executed up to the current state. Since the two-body sintering law affects only the network links, the reference value for a vertex angle is its initial value. The latter need not be the case if a rotational sintering law is introduced, for example, due to asymmetric particle necks. Simultaneous linear equations governing the equilibrium values of the distortions {Auj} and {Atoj} can be obtained. Assume that at the beginning of a sintering time step At the particle network is topologically connected. The two-body sintering transformation is applied to each link for the step At. Quasistatic equilibrium values for {Auj} and {Ato~} are obtained by minimizing equation (2) with {Auj} and {Aeoj} as variables, subject to the condition that the network retain its topological connectedness, i.e. the map of link and vertex interconnections that existed before the time step was executed must be retained. The retention of topological connectedness is equivalent to the requirement that closed polygons in the network remain closed and the sum of angles around each interior vertex remains fixed at 2n. These closure conditions are illustrated in Fig. 6. As outlined in the Appendix, the closure conditions for polygons and vertices can be expressed as a set of linear constraint equations in {Auj} and {Atoj} of the form

)

i

J

where ar bj and cj are constants determined from the

LEU et al.: COMPUTER SIMULATION OF PARTICLE SINTERING

15 0 5 ~ ~ 1

1981

for reconstruction. Using the previously determined map of link and vertex interconnections, additional polygons are fitted into place, as in a jigsaw puzzle, until the network is complete. Because of the closure conditions no gaps can exist between adjacent polygons. The reconstructed network maintains a proper topological connectedness while simultaneously undergoing densification due to the sintering transformation. After reconstruction, newly formed or broken network links can be identified and the linkvertex interconnection map updated for the next time step. The development of specific computer algorithms to search, index, map and reconstruct the network is conceptually straightforward, but nontrivial if optimization is desired. These algorithms involve various pattern recognition and listing procedures that will not be discussed further here.

]2

(a)

Creep mechanisms

(b) Fig. 6. Closure for a network polygon (a) requires that the links ~ defining the polygon edges have a zero vector sum. Closure for an interior vertex (b) requires that sum of the surrounding vertex angles 0j be fixed at 2n.

current network geometry. An appropriate set of such constraint equations Hk applies for each closed polygon and interior vertex in the network, with summations being determined accordingly. The fact that bridging links were excluded from the network topology, and that periodic boundary conditions were used, simplifies the treatment of the constraint equations. Each link is contained in at least one closed polygon, and each vertex appears as an interior vertex. Such simplification is not, however, necessary for the implementation. Equation (2) can be minimized subject to the topological constraints embodied in equation (3) by the method of Lagrange multipliers. Define the auxiliary function 1

2 s

,

1

(4)

where the first two terms are the energy U from equation (2) while the last term is a summation over all constraint equations (3) with {2k} being the set of Lagrange multipliers. As outlined in the Appendix, equation (4) is minimized with {Au,}, {AoJj} and {2k} as independent variables. The equilibrium values of {Aui} and {A%} determined from the minimization of equation (4) permit a proper topologically connected particle network to be reconstructed after any time step. A set of individual transformed, closed polygons are first generated using the {Auj} and {A%} values. Any particular polygon can then be chosen as the starting polygon

Once the elastic distortions {Auj} and {A%} are determined, internal forces and moments acting between particles in the sintering network can be evaluated using the stiffness coefficients. This allows creep mechanisms to be treated within the framework of the simulation. A basic creep effect can be illustrated for the particle pair shown in Fig. I. In addition to shrinkage driven by the sintering mechanism, the force F associated with the topologically necessary elastic distortions provides an additional driving force for changing the link length l by axial creep. In order to treat creep response for the particle pair shown in Fig, 1, a two-body creep model is needed. For given particle-pair geometry and force F, an unconstrained creep displacement Ale = gAt must be evaluated for a time step At. The creep-rate function g = Ale~At is analogous to the sintering-rate function r in equation (1). The function g depends upon the force F. A major advance provided by the current simulation method rests in the fact that it provides a unique, physically-based procedure for evaluating internal forces and moments, hence, creep effects. Creep can be included into the simulation method using superposition. Suppose that AI and AI~ represent the changes in link length I due to sintering and creep, respectively. Then the unconstrained change in link length for the particle pair in Fig. 1 at time step At is 1 ~ 1 - (A/+ A/c). The sign of AIc will depend upon the current force F such that l increases for tensile force and decreases for compressive force. In effect, creep is treated using a superposition of stressfree transformations for sintering and creep. Within the quasistatic approximation, the superposition of kinetic processes such as sintering and creep will be valid in the limit of small time steps. This can be verified explicitly by considering very simple networks and comparing the simulation results to exact solutions of the appropriate differential equations [22].

1982

LEU et al.: COMPUTER SIMULATION OF PARTICLE SINTERING

For the implementation of the simulation code, an empirical, two-body viscous creep model giving the creep law g = t/F

(5)

was used. AIc = g A t = ~lFAt and the sign of Ale is determined from F as described above, r/ is an empirical parameter that scales the relative strength of the creep response with respect to the sintering response given in equation (1). No creep response occurs for r / = 0. Just as in the case of the sintering law, the creep law can be arbitrarily chosen, and it need not be a viscous response. Furthermore, a rotational creep law that depends on the moments M (Fig. 5) could also be included, which would dictate the creep response for vertex angle change. Analogue

Consider the particle network shown in Fig. 2(b) to be made up of moveable pegs (vertices) in the plane joined by axial elastic springs (links). For simplicity, ignore any rotational spring elements for the vertices. The array is initially assumed to be in mechanical equilibrium under the action of the springs, which are in various degrees of tension or compression relative to some specified reference state. Now fix the pegs in the plane so that they cannot be moved. Remove each of the springs of length lj and cut a length AIj from it. This defines the transformation/j:~ lj - A/j for each time step. The values of AIj can represent combined sintering and creep. Reattach each of the transformed springs to their proper pegs, and allow the pegs to move in the plane. The system relaxes into a new configuration, retaining proper topological connectedness, while undergoing a certain amount of densification via the transformation. The relaxed configurations correspond to mechanical equilibrium of the elastic elements, which react instantaneously for quasistatic behavior (no inertial mass), Topological rearrangement thus follows a unique path determined by minimum stored elastic energy. This picture provides a simple analogue for the sintering simulation method. If pegs and springs were arranged in a three-dimensional configuration, a similar analogue applies. From the analogue, it is clear that the simulation method can be extended to three-dimensional particle arrays. In this case, the closure conditions for topological constraint involve network polyhedra and solid angles subtended at each interior vertex. The mathematical formulation is considerably more complicated in detail, but the same underlying physical ideas apply. COMPUTER SIMULATIONS Parameters

Random packings of two-dimensional arrays of particles were generated using the computer algorithm discussed by Hare [12]. A square area with

Fig. 7. Two-dimensional randomly packed particle array (top) and the corresponding network (bottom). The packing density p = 80%. Periodic outer boundaries are used and the unit cell for the array is contained within the dashed lines.

initial dimensions scaled to 1 x I was used as the starting unit cell. Figure 7 shows a typical computergenerated array of particles with periodic boundaries. The unit cell is shown by dashed lines. During sintering particles interpenetrate and the size of the unit cell in the vertical and horizontal directions decreases, although the decrease need not be the same for both directions. After a sintering time step is executed, new values for the edge lengths of the unit cell in the horizontal or vertical directions are determined by calculating the projected path length in that direction of a sequence of links joining the same particle over one period of the array. Because of periodicity, all such paths must give the same edge length, no matter what particles or sequences of links are used.

LEU et al.: COMPUTER SIMULATION OF PARTICLE SINTERING

1983

The (percent) particle packing density p is defined boundary conditions is then the unit cell for the as p = IOONAp/A¢ where N is the number of particles lattice. For ordered arrays, the densification rate in the unit cell, Ap is the particle cross-section area follows directly from the two-body sintering law and A¢ is the current unit cell area. Ap is fixed by the because the unit cell shrinks uniformly, without topoinitial particle size so that an increase in p due to logical constraint. Furthermore, the symmetry reparticle interpenetration is embodied in the reduction quires that the elastic distortions be identically zero. of Ac. No correction is made for particle overlap The simulation code produced these specific solutions geometry. The starting arrays for simulations con- when applied to ordered arrays of particles. tained about 50 particles with p = 75-80% (the maximum p value is 91% for planar hexagonal arrays). Simulation results It was necessary to keep p in this range in order to The simulations were made using double-precision eliminate bridging links (Fig. 3). Loosely packed particle arrays in the range p < 75% are thus ex- Fortran on a VAX 11/750 computer with 3 Mb of cluded from the initial implementation of the simu- memory. Details on code development are available elsewhere [22,23]. The Sparspak [24] subroutine lation code. The stiffness coefficients k and 13 were taken to be package was used to solve the linear minimization constants for the simulation runs. The numerical equations. A typical sintering time step took several value of k was normalized to unity so that specifying minutes to execute with about 50 particles in the unit values of fl is equivalent to adjusting the ratio fl/k. cell. Apart from the Sparspak routines, the code was The numerical values of fl used for the simulations used in first-generation form and no concerted effort runs were 10 -6, 10 -3, 1 and 103. As is shown in the was made to optimize the various algorithms. SubAppendix, a normalization of equation (2) implies sequent experience indicates that significant imthat links and vertex angles in the network will have provements in execution time can be achieved with a comparable elastic "flexibility" when fl is about 10 -2. moderate amount of optimization. Figure 8 shows a typical result for the topological The strength of the creep response was determined by the scaling parameter ~/ in equation (5). Simu- rearrangements that occur in a unit cell during a lations were run with r/values ranging from 0 to 50. simulation run. Densification is manifest by the decrease in the unit cell size. As can be seen in the figure, Verification qf the simulation code major particle rearrangements do not occur, the Detailed verification of the numerical results pro- densification being reasonably uniform and isotropic. duced by the simulation code is a difficult task However, more subtle particle rearrangements occur because there are very few exact solutions against that cannot be easily detected in Fig. 8. The rather which it can be checked. Nevertheless, several im- large particle rearrangements observed in various portant criteria arise when implementing and veri- experiments [2-5], and in the earlier computer simulation studies of Ross et al. [9], were not observed fying the simulation code. An efficient choice of the time step At for the here. This is attributed to the fact that p for the simulation runs is important. This can only be done starting arrays used here was kept in the range of by a trial and error procedure. Experience with the 75-80%. The absence of bridging links (chains) in the version of the simulation code implemented here, starting arrays at higher packing densities no doubt with the choices of sintering, creep and stiffness inhibits the propensity for particle rearrangement, as parameters already mentioned, resulted in a value can be inferred from earlier results [2-5, 9]. SomeAt = 0.02. Using this At, relatively rapid densification what more particle rearrangement could be achieved could be achieved while the distribution of forces and when the relative angle flexibility was high moments in the links and vertices of the network ([3 = 10-6), however, even in this case extensive rearchanged in a smooth, continuous fashion. If At is too rangements did not occur. Figure 9 shows the changes in the total stored large, numerical inaccuracies or oscillations may ocelastic energy U and in the total number of links in cur. After the minimization step, the numerical closure the unit cell during a series of sintering time steps for errors for polygons and vertex angles were typically five different starting arrays. The iteration number the order of 10 -~-'. This approaches the overall nu- identifies the time step. For these runs, [t = 1, r / = 0 merical accuracy for the computations, using double and the nominal starting density ,o = 80%. U inprecision. The edge lengths for the unit cell remained creases as sintering progresses, which was a feature constant to at least the order of 10 6 for any given common to the simulation runs for all starting conditime step when averaged over many different particles tions. An increase in the number of links during and link sequences. This level of numerical accuracy sintering implies that new contacts have been formed between particle pairs. A comparison of the curves in is considered acceptable. The only "benchmark" solution for the sintering Fig. 9 shows that the rather abrupt increases in U can problem in large-scale particle arrays is that for be correlated with the formation of new contacts. The perfectly ordered arrays, i.e. square or hexagonal densification behavior for the different starting arrays particle lattices. The network unit cell for periodic in Fig. 9 was very similar.

1984

LEU et al.:

COMPUTER SIMULATION OF PARTICLE SINTERING

10 8 >. O = W Z m

6

= W

4

I,r

2 10

5

15

20

iS

20

ITERATION

80"

"7-

60. Z

50

o

S

lb ITERATION

Fig. 9. Stored elastic energy U (arbitrary units) vs iteration number for five different starting arrays with N = 50, fl = l, t / = 0 and a nominal packing density of p = 80% (top). Total number of links in the unit cell for the same arrays (bottom).

Fig. 8. Starting array (top) for a simulation run where N = 46, p = 73%, fl = l and rl = 0. Final array (bottom) with p = 90% after a sequence of sintefing time steps. The array unit cell is shown by the dashed lines. Figure 10 shows the change in U for a series of sintering steps for different values of ft. The same starting array was used for all of these simulation runs. When the value of fl increases from l0 -3 to l, there is a significant increase in the magnitude of U. As was noted earlier, the links and angles in the particle network have comparable flexibility when fl is about l0 -2. For smaller fl values the angles are more flexible while for larger values the links are more flexible. This suggests that when fl > l, elastic distortions are accommodated primarily by the more flexible links and the values of U approach a limit determined by the fixed link stiffness k = 1. When fl < l, elastic distortions are accommodated to an increasing extent by the more flexible angles and the value of U decreases with ft. A different representation of the stored elastic energy in the particle array is given in Fig. I I. A local energy value can be assigned to each network vertex by adding the energy due to each of the surrounding vertex angles plus half the energy due to each of the links joining that vertex. At any point in the sintering history, a m a p can be drawn of the energy distribu-

tion or fluctuation in the particle array. Figure 11 shows typical maps of this kind at the 10th iteration for the simulations in Fig. 10 for values of fl = 10 -6 and 103. The circles drawn at each vertex in the network are scaled according to the lcoal energy value for that vertex. For clarity, the m a x i m u m circle size in the diagram is held fixed so that absolute comparisons between diagrams should not be made. Figures 10 and 11 shows that changing fl not only changes the magnitude of U, but also its distribution. As sintering progresses the distribution undergoes smooth change. N o specific correlations could be drawn between energy distribution and local particle

.=.

4 3

¢~

2

¸

0

S 5

10

15

20

ITERATION

Fig. 10 Stored elastic energy U (arbitrary units) vs iteration number for a starting array with N =50, ~ = 0 and p = 80%. The same starting array was used for different simulations using the values of fl indicated.

LEU et al.: . , 4 . . . .

•

.

COMPUTER SIMULATION OF PARTICLE SINTERING

•

..

,...::" '

• .:•

'

•

,

••

: .... i . Q

•'

.0

'

:

'

'

• ' .... : .•..~...

•

.,•'.

,

B

,

,

'.c

' • • ••

",="

,

.=."

. . . . :,-:••:•

.... @ •

"a

,

' ,.

,

.':,

:.

~ , '

. -•°,

,

"

:':

,

,

•••

•......

."

'• ..

"

,..'J.

-

......

o

'.¢

.,4

.

, • .

,

-

•' •

,

, .,,~ ~0• -• . -

DISCUSSION .

, * • .':

..••~...

ii!!ii !! ,

,

•..

"

. tp,

..

.-.

•

Fig. 11. Energy fluctuation maps as described in the text for 10-6 (top) and fl = 103 (bottom)• For clarity, the maximum size of the energy circle is the same in all diagrams even though the maximum vertex energy values may differ• fl =

geometry. In general, it appeared that more regularly packed regions of the array can be associated with lower energy• Figure 12 shows the change in density p during a series of sintering time steps for different values of the creep parameter ~/. The same starting array with fl = I was used for all of the simulation runs. The slope of the lines represents the densification rate. Relatively constant densification rates were observed in all simulations, except for certain starting arrays at lower initial densities where variations could be observed. The creation of new particle contacts appeared to play an important role in this regard. The line for ~/= 0 in Fig. 12 is representative of the densification behavior for the different arrays corresponding to Fig. 9. AM.

361~H

With increasing values of r/ in Fig. 12, the densification rate decreases such that at r / = 50 there is a significant retardation in rate When q = 50, the average changes in link length due to creep become comparable to the changes due to sintering If r/ is made too large, the changes in link length for a time step At = 002 become excessive and there can be numerical instabilities arising during the simulation run. Reductions in At permit larger values of r/to be used, however, the results shown in Fig 12 are sufficient to illustrate the basic creep interaction Energy distribution maps of the kind shown in Fig 11 showed that when creep is present, the distribution of stored elastic energy becomes more uniform, i e creep tends to smooth out the fluctuation amplitude. All of the creep response shown in Fig. 12 was due to axial creep Further synergism due to rotational creep effects can certainly be anticipated.

:.,...~,...

- .~

.o..

1985

One of the primary goals of the current investigation was the development of a physicallybased simulation method for early stage sintcring along With an initial implementation of a computer code to demonstrate the viability of the approach• The simulation results obtained were parametric in nature and were based on very simple empirical models• The results provide proof-in-principle for the effectiveness of the method; they are not intended as a study of kinetics and mechanisms for a particular material system. The computer simulation method developed here provides an effective framework for treating topological constraints and stress-related effects. The method is quite general in the sense that arbitrary sintering or creep laws can be introduced. Both axial and rotational effects can be included• Relevant elastic properties are introduced through stiffness coefficients that uniquely characterize the elastic properties of the particle network. Creep effects can be included via a separate creep model, as done here, or they may be included via a stress-dependent,

.92 q=lO rl=O

Q

.80

0

lo

2;

3o

ITERATION

Fig. 12. Density p vs iteration number for a string array with N = 50 and fl = 1• The same starting array was used for different simulations using the valucs o f r/indicated.

1986

LEU eta/.:

COMPUTER SIMULATION OF PARTICLE SINTERING

two-body sintering model. Several models of the size distribution. As far as the current authors are latter kind have appeared in the literature [1 l, 25]. aware, the simulation result provides the first direct Finally, with suitable modification, the simulation verification that such intrinsic sintering-creep synmethod can be extended to three-dimensional particle ergism exists. The significance of the result lies in the fact that the actual densification kinetics for real arrays. Extensive particle rearrangement was not observed material systems can be a complex interplay between during the sintering simulation runs because the the basic sintering mechanisms and creep mechsimulations were made with starting arrays having anisms driven by topologically-induced stress. The mechanistic origin of densification retardation relatively high packing densities. Attempts to start with low density arrays invariably resulted in link due to creep is difficult to analyze. A detailed exambreakage during sintering and subsequent production ination of the simulation results showed that when of bridging links and chains (Fig. 3). While this creep is present, there is a complicated distribution of behavior could not be treated with the current imple- internal forces with some network links remaining mentation of the computer code, it does point to under tension and some under compression. Thus, on significant and important rearrangement effects at a local scale, the axial creep response can either aid lower densities. Nevertheless, more subtle particle or retard densification. The net effect must then be a rearrangements did occur in the denser arrays and statistical one, where retardation dominates when these were dependent on the value of the stiffness appropriately averaged over the network. An interparameter ft. For large values of fl, densification pretation in terms of a long-range hydrostatic stress occurred by an almost homogeneous shrinkage of the component [20] is not easily adapted to the intrinsic network polygons with their angles fixed. For smaller effect observed here. In summary, the viability of a new physically-based values of /~, the vertex angles become relatively flexible and more rearrangement occurred due to computer simulation method for early stage sintering polygon shearing. This suggests that the relative has been demonstrated. One of the key features of the angular stiffness fl plays an important role for particle method is the recognition that topological constraints rearrangement in loosely packed arrays. Recent sta- in large-scale random particle arrays can lead to tistical studies of particle rearrangements for two- elastic distortion and associated stress. The stiffness dimensional arrays of sintering particles by Weiser parameters that describe the elastic response of the and DeJonghe [26] shows that vertex angle change is particle array play an important role in rearimportant for sintering. Some simulations were made rangement effects. Furthermore, the elastic stresses here using small clusters of 5-10 particles with free can drive a creep response that results in an intrinsic (non-periodic) outer boundaries. These small clus- sintering-creep interaction and a retardation of the ters, while lacking bridging links or chains, are still densification kinetics. relatively unconstrained at higher packing densities. Acknowledgements--The authors gratefully acknowledge Noticeable particle rearrangement during sintering the Department of Materials Science and Engineering, could be observed in these clusters, especially at North Carolina State University, for support of this resmaller values of fl, i.e. for high vertex angle search work as well as for providing the necessary computer time. This work was performed in partial fulfillment of the flexibility. Ph.D. requirements for H. J. Leu. The densification rate was not significantly influenced by changes in the network topology reREFERENCES suiting from different starting arrays. This again 1. W. S. Coblenz, J. M. Dynys, R. M. Cannon and R. L. appears to be a manifestion of the higher packing Coble, Sintering Processes, Mater. Sci. Res. (edited by densities used since Ross et al., observed certain G. C. Kuczynski), Vol. 13, p. 141. Plenum Press, New topological influence effects for loosely packed parYork (1979). ticle arrays [9]. 2. H. E. Exner, G. Petzow and P. Wellner, Sintering and Related Phenomena (edited by G. C. Kuczynski), p. 351. One of the important results obtained in the initial Plenum Press, New York (1973). simulation studies deals with sintering-creep inter3. H. E. Exner and G. Petzow, Sintering and Catalysis, action. When creep occurs during sintering, the Mater. Sci. Res. (edited by G. C. Kuczynski), Vol. I0, densification rates noticeably decrease (Fig. 12). A p. 279. Plenum Press, New York (1975). related effect has been observed experimentally using 4. G. Petzow and H. E. Exner, Z. Metallk. 67, 611 (1976). 5. H. E. Exner and P. Bross, Acta metalL 27, 1007 (1979). composite mixtures of powder particles. When a 6. A. Nayala, L. Mansur and J. White, Powder Metall. 6, small volume fraction of large non-sinterable par108 (1963). ticles are introduced into a matrix of sinterable 7. R. L. Eadie, W. A. Miller and G. C. Weatherly, Scripta particles, the densification rates are significantly remetall. 8, 755 (1974). 8. H. E. Exner, Rev. Powder Met. Phys. Ceram. I, 1 (1979). duced [16]. The non-sinterable particles produce 9. J. W. Ross, W. A. Miller and G. C. Weatherly, Aeta internal stresses due to topological constraint, thus metall. 30, 203 (1982). driving creep processes that act to retard densification 10. P. Bross and H. E. Exner, Acta metall. 27, 1013 (1979). [18-20]. However, the retardation observed in the I 1. J. W. Ross, W. A. Miller and G. C. Weatherly, J. appl. current simulation study is an intrinsic effect that does Phys. 52, 3884 (1981). not depend upon differences in particle properties or 12. T. M. Hare, Sintering Processes, Mater. SCi. Res. (edited

LEU et al.:

13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 26.

COMPUTER SIMULATION OF PARTICLE SINTERING

by G. C. Kuczynski), Vol. 13, p. 77. Plenum Press, New York (1979). N. Ramakrishnan, T. Balakrishna Bhat and V. S. Arunachalam, Acta metall 32, 357 (1984). J. W. Ross, Ph.D. thesis, Univ. of Toronto (1980). F. F. Lange, J. Mater. Res. 2, 59 (1987). L. C. DeJonghe, M. N. Rahaman and C. H. Hsueh, Acta metall. 34, 1467 (1986). C. H. Hseuth, A. G. Evans, R. M. Cannon and R. J. Brook, Acta metall. 34, 927 (1986). C. H. Hseuh, A. G. Evans and R, M. McMeeking, J. Am. Ceram. Soc. 69, C64 (1986). C. H. Hseuh, Acta metall. 34, 1497 (1986). R. Raj and R. K. Bordia, Acta metall. 32, 1003 (1984). D. J. Bacon, D. M. Barnett and R. O. Scattergood, Prog. Mater. Sci. 23, 51 (1978). H. J. Leu, Ph.D. thesis, North Carolina State Univ. (1987). H. J. Leu, Program Listing Supplement, Ph.D. thesis, North Carolina State Univ. (1987). A. George, J. Liu and E. Ng, User Guide for Sparspak. Waterloo Sparse Linear Equations Package. Dept. Computer Sci., Univ. of Waterloo, Ontario, (1980). C. J. Cosgrove, J. A. Strozier and L. L. Seigle, J. appl. Phys. 47, 1258 (1976). M. W. Weiser and L. C. DeJonghe, J. Am. Ceram. Soc. 69, 822 (1986).

APPENDIX

H~= ~ ( - I) ~+~(l~ + aU~)COS0~

? Au,

?W ?A~o,

- - = 0

?W = 0. (A4) 72, A system of simultaneous linear equations results from these derivatives. Again, without producing the detailed analysis, the system can be represented in block matrix form by D A X

,A,,

where D is a diagonal matrix containing the stiffness coefficients k i and fit. A and its transpose A r and • are matrices containing constants determined by the current network geometry. ~ - [ A u t . . . . . A(o, . . . . ] is a column vector containing the elastic distortion variables and 2 = [ 2 t .... ] is a column vector containing the Lagrange multipliers. 0 denotes a zero matrix or vector. A solution to the minimization equations is obtained by expanding the block matrix system in equation (A5) into the two sets of matrix equations

(--1)J÷'l, sin0tAwk=0

A T D-i A2 = - B

(A7)

Further consideration of this system will show that the square matrix ATD-~A is sparse. Linear equation solvers like the Sparspak routines [24] are very well suited to this system. Once the solution for 2 is obtained, the desired solution for X is obtained by back substituting the result into equation (A6) (AS)

(AI) Network Flexibility

H 2 = ~ ( - 1)3+ I(l) + Auj) sin 0j j=l n~l

~

(A6)

Since D is diagonal, the inverse D -~ is readily obtained. Equations (A6) can then be combined to obtain a determinate system of equations for the Lagrange multipliers 2

X = -D-IA2

n

~

k=ly=k+t

+Z

- - = 0

A TX = B.

The network shown in Fig. 2 will be properly connected if each minimal polygon (a polygon that cannot be reduced to another polygon by eliminating links) remains closed after a time step At and the sum of the angles about each interior vertex (a vertex completely enclosed by minimal polygons) remains fixed. Without providing the straightforward but lengthy derivations, the closure conditions for the n-sided polygon shown in Fig. 6(a) are given by (horizontal and vertical vector components)

n~l

The minimization equations are obtained from the derivatives ?W

D)~ + A2 = 0

Minimization Equations

+~

1987

( - l ) J + ~ l j c ° s 0 , A(°'=0"

(A2)

k~ l j=k+l

For the polygon shown in Fig. 6(a), n = 5 and the labeling scheme is indicated. The link lengths Ij and angles 0j are the stress-free reference values for the current time step. In deriving equations (AI) and (A2), it was assumed that the changes Acoj are small. Since the sum of angles around a vertex must remain fixed at 2n, the closure condition for an m-angle vertex is m

H 3= ~ A w / = 0

(A3)

/=1

For the vertex shown in Fig. 6b, m = 4 and the labeling scheme is indicated. Equations (AI)-(A3) have the general form shown by equation (3) in the text where the summations and constants aj, bj and cj are to be determined for each polygon and vertex in the network. The constraint equations are linear in the variables AIj and A~o). Equation (2) in the text can therefore be minimized subject to the constraints in equation (3) by the method of Lagrange multipliers. As was indicated in the text, this lollows by defining the auxiliary function W in equation (4) where {Aui}, {Awi} and {2,} are the independent variables.

An estimate of the relative effects of the stiffness coefficients k and f on network flexibility can be obtained using the stored elastic energy U given by equation (2). Let ky = k and pj = f for all values ofj. Further, let L and (9 be average values for the distribution of link lengths {//} and vertex angles {0r }, respectively. The summations in equation (2) can be written in the form of sums over nondimensional strains as follows

As a first approximation, it will be assumed that links and angles make about equal contributions to U when the multiplying factor for the second sum is unity. Therefore, comparable contributions to U as a result of network flexibility obtain when B

L: = ~.

(B2)

A numerical value for the critical ratio #/k depends on L and (9. There will be the order of 6 vertex angles surrounding each vertex (hexagonal packing) and (9 ~ 2rr/6 ~ I. Likewise, L 2= Ac/N so that for the conditions used here with about 50 particles in a 1 x I starting unit cell, A~ = 1, N = 50 and L: = 0.02. Thus, the critical value f / k ~ 10 -2 as noted in the text.