3-DIMENSIONAL STRUCTURE-IMITATION MODEL OF EVOLUTION OF MICROSTRUCTURE OF POWDER BODY DURING SINTERING

This paper describes a 3D structure-imitation computer model ofevolution ofthe powder compactduringsinteringand recrystallizationwithoutnucleation. Atthe initial stages ofthe evolution processes (sintering until a mosaic structure ofboundaries is formed) the model particles are spheres, and two-particle interaction laws control their evolution. During sintering the degree ofmutual penetration ofthe particles increases, regions where spherical particles are wholly facetted by contacts with neighboring particles are formed and grow. These particles are described using the formalism ofVoronoi radical polyhedra, and grain growth laws govern their evolution. The model predicts the time dependencies of the following structure parameters of the polyhedra: average polyhedron size and dispersion, total surface ofthe facets of the polyhedra and total length ofthe edges of the polyhedra.


INTRODUCTION
Analyzing the physical phenomena that occur during sintering in the powder body, it turns out that the evolution of the microstructure plays the main role (Geguzin, 1984;Skorokhod and Solonin, 1984; Kadushnikov and Beketov, 1989;Ross et al., 1982).The existing models of microstructural changes can be reduced to two main types, which are based on pair-interactions and on the evolution of the topological structure.For the first case an advanced apparatus to describe the kinetics of the process through a phenomenological equation of the sintering theory has been developed (Geguzin, 1984).Within the framework of this approach the elementary sintering models of the system consisting of equisized spheres are investigated.In the second case model approximations of the microstructure of the material by sets of polyhedra are constructed.In agreement with some main principles (maximum entropy, equilibrium of dihedral angles, etc.) imitation models of grain structure evolution are investigated.In all cases, simulations are performed for separate stages of the entire sintering process.
The authors have undertaken the approach to construct a uniform concept ofsimulation ofthe microstructure evolution during all stages of the sintering process (Kadushnikov and Beketov, 1989;Kadushnikov and Alievsky, 1989;Kadushnikov et al., 1991 a,b,c;1998).The sintering process is described by 4 main stages.Each stage is characterized by a different geometrical model of the microstructure (Table I).
The first 2 stages correspond to the initial stages of the sintering process.A geometrical model of the microstructure is a set of spheres.
The 3rd stage corresponds to the late stage of the evolution of the sintered bodies, when the degree of mutual penetration of particles is high enough and the quantity of completely closed particles becomes significant.The set of spheres is still used for mathematical simulation of the sintering process.However, the skeletal polyhedral and triangula- tion models (corresponding to the configurations of spheres) are being Initial stage of the sintering: mutual penetration of the spheres is low Mutual penetration of the spheres is high The sintered bodies occupy the whole space and are described by polyhedra built during the sintering process.Here skeletalpolyhedralmodelis based on the 3D Voronoi diagram built on the given set of spheres, skeletal triangulation model is based on the 3D Delaunay triangulation corre- sponding to this Voronoi diagram (Preparata and Shamos, 1985).These models describe the real physical processes better than the set of spheres at the stage 3. Terms corresponding to the development of grain structures in the phase of normal grain growth are added to equations of the dynamics of the sintering spheres.
The 4th stage corresponds to the final phase of the sintering process.The set of spheres is not used now.Instead, the skeletal models built on the last step of the stage 3 become the base geometrical models.Fol- lowing evolution of the system is determined by the dynamics of the skeletal models.This stage corresponds to recrystallization processes.
In this paper, we give a brief description of the algorithms used to simulate all stages of the sintering process.To build a close packing of spheres with disperse sizes within a bunker with flat walls we chose an algorithm named "rolling algorithm".The physical idea of this algorithm is formulated as follows.The spheres are generated according to the chosen law of size distribution and are dropped to the bunker either from one point, or from randomly chosen positions.As soon as the dropped sphere encounters an obstacle the bunker wall or an already packed sphere it sticks to it (without impact) and begins to slide on its surface in the direction of the minimum of the potential energy to the following obstacle.(This direction is a projection of the free fall direction on the surface of the obstacle.)The movement of the sphere stops, obviously, at a point of intersection of 3 surfaces (3 spheres, 2 spheres and plane and etc.) or on a surface situated perpendicular to the direction of the free fall of the sphere (for example, on the bottom of the bunker).

DESCRIPTION OF THE MODEL
It should be noticed that the movement of the current packed sphere is completely described by the trajectory of its center.In this case, the sliding of the sphere with center O and radius R on the surface of another sphere with center O1 and radius R1 corresponds to sliding of the center O on the surface of a sphere with center O lJand increased radius R1 + R. Similarly, the sliding of the sphere with center O and radius R on the plane ax=b, (1) a (ax, ay, az), a > 0, corresponds to sliding of the center O on a shifted plane ax b + laiR.
(2) Thus, it is possible to formulate the following algorithm of movement of the next packed sphere: (1) Choose radius R of the packed sphere and an initial position of its center O. Initialize the set of restrictions: no restrictions.
(2) Analyze the current set of restrictions: the system of 0, 1, 2 or 3 inequalities describing positions in the space which are disabled for the center of the packed sphere.Choose the new direction of movement of the sphere center.The movement of the center along this direction must satisfy the restrictions, this direction must form the minimal allowable angle with the direction of the free fall of the sphere (vector (0, 0, 1)), and this angle must be less than r/2.In absence of allowable directions, i.e. when all angles are more than r/2, stop and choose the next sphere for package.
(3) Determine the next segment of the movement trajectory of the center.It is a direct line or an arc segment.Perhaps, remove one,of the restrictions.
(4) Calculate points of intersection of the current segment of the tra- jectory with all packed spheres (with radii increased by R) or bunker walls (shifted by a normal vector with length R).Select the nearest to the beginning of the segment.Add one new restriction: all points inside the sphere or behind the wall (which intersects the current segment at the nearest position) are disabled for the center of the packed sphere.Return to step 2.
It is easy to see that the described algorithm depends quadratically on the number of packed spheres.To convert it to a nearly-linear algori- thm, the special technique of grid hashing was used.This technique will be outlined in the following.
Let us consider a cubic grid with the cell size comparable with sizes ofspheres.For each cell, let us build the list ofpacked spheres the centers of which lie in this cell.Then analysis of spheres disposed in several grid cells around the trajectory allows to find all intersections of the trajectory segment and the packed spheres.So, when the parameters of the grid are chosen well, the algorithm becomes almost linear with a multiplier depending on the parameters of the whole packing of spheres (such as dispersion.ofradii of the spheres, etc.).
Stage 2: Evolution of the Packing of Spheres The constructed close packing of spheres represents the initial data for the second stage simulation of the evolution of the spheres during sintering.For deriving the equations of movement and dynamics of the sizes of the sintered body, the main incentive was to satisfy the equations of kinetics of pair sintering for each pair of contacted spheres in the system.
Let us consider sintering of a system consisting of N spheres: Si (Ri, rg, vg), i-l, 2,..., N; Ri radii of the spheres, ri coordinates of their centers, vspeed of movement.
To describe the free movement of the spheres, the well-known equa- tions of Newton are used: i, miigi-" Fi. (3) Here mi -rpR, p density of the material of the spheres, Fi sum of external forces, acting on each sphere, such as the force of gravity (m/g) or force of liquid friction (-v vi[S[i, [S[i area of free surface of ith sphere).
At existence or occurrence of contact: Iryl-Iri-rjl < Ri + Rj, (4) it is necessary to take into account the interaction of the contacted spheres, resulting in sintering.Let us record the equations of the kinetics of pair sintering for all pairs of incident (contacting) spheres in the following form: v g =f(lrl)r;, i= 1,2,...,N; j---, i; (5) here and belowj means that the spheres andj are contacted see (4).
And let us construct the deviation functional: Q Z Ivo'-f(Iril)ryl " (6) i,j:j---* Let us modify the trajectory of movement of the corresponding point for an ensemble of spheres in the phase space so that it comes nearer to a minimum of the deviation functional Q.To do it, let us add the anti- gradient Q with a weight multiplier a in the right parts of Eq. ( 3), where the multiplier a determines the degree of influence of sintering on movement of the system of spheres (or vice versa): The conservation of weight means that during movement and sin- tering the volume V of the geometrical union of spheres remains constant: OV. OV ) + o. (81 Let us balance the change in volume V of the union of spheres con- nected with the movements of some sphere Si with the change in the radius Ri of this sphere: OVi + OV Or---./ 0 (9) for each i.Thus, we obtain the following equations ofdynamics ofthe radii ofthe spheres: ke-Ielze/ISel, (10) Here ISel ff , ds (OV)/(ORi) area of the free surface qi of the sphere S, .'-.i ff,,(Mll) as can be interpreted as an oriented area of the projection of occupied surface i of the sphere Si on the plane per- pendicular to its velocity (n is the unit normal vector to the surface).
The integration of the systems of Eqs. ( 7), ( 10) is carried out by the Euler method with fixed time step (Hockney and Eastwood, 1981).Areas ISil, and .=.i, are calculated numerically.The packing of 1000 spheres with distributions of radii N(1.0, 0.3) in a cubic bunker and sections of the model during the evolution process are shown in Fig. 1.

About Skeletal Models
To understand the algorithm used to simulate the sintering process at the stages 3 and 4, it is necessary to describe the concept of skeletal models and some connected terms.
We consider 2 kinds of skeletal models based on Voronoi polyhedra and 319 Delaunay triangulation.These geometrical concepts are some generalization of the classic concepts described, for example, in Preparata and Shamos, 1985.These models allow to transit naturally from the set of spheres to polyhedral description of the system.
The basic concepts of our skeletal models will be outlined in the following.
Let there are the set ofN spheres: Si (Ri, ri), 1, 2,..., N; Ri radii of the spheres, rradius-vectors of their centers.We define the Voronoi polyhedron P as a set of points {r R 3" Ir ril R/2 _< [r r[ Rj 2. for all j = i} (11) (it is equivalent to the classic definition when all R 0).Voronoi poly- hedra occupy the whole space R3.Adjacent pairs of polyhedra Pi, Pj are separated from each other by radicalplanes ofpairs ofthe corresponding spheres Si, Sj {r R 3" [r-ri[ 2 R/ [r rl z Rj2.} (12) (when the spheres Sg and S are intersected, their radical plane contains the intersection circumference).So, we also call the Voronoi poly- hedra "radicalpolyhedra".We also will consider radical axis of 3 spheres s,, sj, s: R} (13) and radical center of4 spheres S, S, S, SI: the point r determined by the condition 2. (14) Delaunay triangulation is defined as the partition of the convex hull of the set of centers of all spheres Si into tetrahedra, such as: vertices of the tetrahedra are all centers r; there is an edge of some tetrahedra connecting the vertices r and rj if (and only if) the radical polyhedra Pi, P are adjacent (have a common facet); there is a facet of some tetrahedra with the vertices r, r and rk if (and only if) the radical polyhedra Pi, Pj, Pk are adjacent (have a common edge); there is a tetrahedron with the vertices ri, rj, rk and rtif(and only if) the radical polyhedra Pi, P, Pk, Pt are adjacent (have a common vertex).
Our skeletal triangulation model is obtained from the Delaunay tri- angulation by removal of edges corresponding to not contacted spheres ([ri-rj[ 2 > (R n t-Rj)2) and removal of all facets and tetrahedra contain- ing these edges.Similarly, our skeletalpolyhedral model is obtained from the Voronoi polyhedra by removal of the corresponding facets, edges and vertices.
The duality of the triangulation and polyhedral skeletal models is convenient for description of grains by the radical polyhedra and pores by tetrahedra of the triangulation.
It is assumed that the radical plane of each pair of spheres separates their centers.This property is kept during the evolution of the sintering of the spheres.
Both skeletal triangulation model and skeletal polyhedral model can be trivially built ifwe have already built the Delaunay triangulation.So, we give a brief description of an algorithm of building Delaunay tri- angulation in the following.
This algorithm is iterative.Let some tetrahedra of the triangulation be already built.The algorithm chooses any external facet of some tetrahedron Tm ("external" means that it is not a facet of any other tetrahedron) so that ff does not belong to the convex hull, and builds a new tetrahedron T of the triangulation on the facet ff (the second tetra- hedron incident with this facet).Vertices of the facet ff are centers of 3 spheres Si, Sj and Sk.Building the tetrahedron T means selecting the fourth sphere St the center of which will be the fourth vertex of T. The algorithm selects the sphere St from all spheres lying."above" the plane of the facet (where "above" means "from the other side than Tm") so that the position of a radical center of 4 spheres Si, Sj, Sk, St is the most "lower".Formally, argmin zrp";:-c),:--c where c is an intersection point of the plane of the facet I, and a radical axis of spheres Si, S, Sk ("radical center of the facet "), nc the unit normal vector to the facet ff directed "upward", wp Ic-rpl 2-R 2, Wc -le r,I R: (-le ryl 2 Rj:.--le-rl 2 R).
After choosing the fourth vertex of the new tetrahedron T, new edges and facets are added to the triangulation.The process of building is finalized when all external facets oftetrahedra belong to the convex hull.
This algorithm, as well as the "rolling algorithm" from stage 1, can be made almost linear by means of the grid.
Stage 3: Simulation Using "Skeletal" Models In the third stage ofthe simulation, the dynamics ofthe skeletal models is practically completely determined by movement and sintering of the spheres.Thus, it is necessary to continue the integration of the system of equations of dynamics of the sintering spheres, building the skeletal models repeatedly and calculating the target characteristics ofthe model.
The influence of effects of normal grain gi'owth is reflected only in the equations of sphere size dynamics for spheres with completely closed (by neighbors) surface (not influencing on conservation of mass) the equation of dynamics ofnormal grain growth is used instead of Eq. ( 10): The plots in Fig. 2 show the time dependencies of the following characteristics of skeletal models: average grain size i.e. the diameter of a sphere, the volume of which equals the volume of the polyhedron; dispersion of polyhedron sizes; total area of the sides of the polyhedra; total length of the edges of the polyhedra.

Stage 4: Simulation of Recrystallization and Normal Grain Growth
The fourth stage is characterized by a change of the basic geometrical model particles become polyhedra without spherical surface.The evolution of the model is determined by minimization of the redundant free energy of system: A aS + eL + "y V, (17) where a, e, 7 denote the specific redundant free energies of facets, edges and vertices of the polyhedra and S, L, V denote the total surface of facets, edge lengths and vertices number, respectively.The model of Voronoi polyhedra, was chosen as a model of particles.As independent variables used in the simulation equations we choose the coordinates and radii of the spheres determining the Voronoi diagram, i.e.A A(ri, Ri).
Unlike the stages 2 and 3, here the set of spheres is used only in mathe- matical formalism and is not essential.The equations of movement of the spheres can now be written as follows: (18) (, b factors ofrnobility).While calculating the right parts we consider all the edges of triangulation, that contain vertices r, and then for each edge all the facets, that contain this edge.Each facet belongs to the pair of tetrahedra, whose radical centers are vertices of the radical polyhedra, forming the edge ofpolyhedron.Let us designate them by c+ and c_, let c) denote the center of the facet of the polyhedron, that corresponds tojth edge.The gradient of A is written as follows: OA Ori OA ORi "i -t" " --C e r fe Ori where Sf 1 / 2 (Ic+ cil2le_ cil 2 (c+ ci, c_ ci)2) 1/2. (20) The summation in ( 19) is performed over all the edges e containing the vertex ri and over all the facetsfcontaining the edges e.The equations of movement for spheres are integrated by Euler's method, with recalculation of skeletal models during each time step.

CONCLUSION
This paper describes the complete concept of simulation of the evolution of the microstructure of the powder compact during all stages of the sintering process.The four stages of the process are considered, the appropriate geometrical model is suggested for each stage.It was shown that the evolution of the system can be described by the equations of movement of the centers and radii of spheres at all the stages of the process.The original algorithms for building the packing of spheres and the skeletal triangulation are presented.The described algorithms for the first three stages are implemented in a computer program.This program is tested with models containing up to 1,000,000 spheres for packing and up to 30,000 spheres for triangulation.We used the personal computer with the Pentium 200 CPU and 64Mb RAM.

Stage 1 :
Building the Close Packing of Spheres FIGUREView of initial spheres package (a), and sections of model during evolu- tion (b) 50, (c)80, and (d)    110 time steps.

FIGURE 2
FIGURE 2 Time dependencies of average polyhedron size (a), polyhedrons' sizes dispersion (b), total surface of polyhedrons facets (c) and total length of polyhedrons' edges (d).

TABLE
Stages of sintering and geometrical models of the microstructure for Stage 2