Direct 3D Simulation of Powder Sintering by Surface and Volume Diffusion

Within the general context of solid-state sintering process, this work presents a numericalmodelling approach, at the grain scale, of ceramic grain packing consolidation. Typically, the sinteringprocess triggers several matter diffusion routes that are thermally activated: surface, grain boundaryand volume diffusions. Including this physics into a high-performance computing framework wouldpermit to investigate and to track the changes occurring into a granular packing during sintering. Inperforming this kind of simulations, one will face several challenges: the strong topological changesappear during sintering simulation at the grains scale, the evolution of the structure is mainly drivenby the surface tension phenomena through the Laplace's law, and the mechanical properties of thegrains could, possibly, be different. The proposed numerical simulations are carried out within anEulerian Finite Element framework and the Level-Set method is used to cope with changes in themicrostructure. The results obtained with this numerical strategy are compared with success to theusual geometrical models.


Introduction
Sintering has became one of the most important manufacturing processes. Sintering is used for the fabrications of high performances materials and parts in a wide -and still growing-range of domains. Since the 50's a large amount of work has been done concerning the theoretical and experimental study of this process. However the understanding of the underlying physical phenomena is still a field of active research. Diffusion phenomena is considered to be responsible of the microstructural evolution that takes place during sintering [1]. After the powder has been shaped, under the heat action different diffusion mechanisms are thermally activated. During the first stages of the process the necks between the particles is created. Necks continue to develop allowing the open porosity to became closed to finally end up by disappear at the end of the sintering process. Even if it is possible to distinguish at least six different diffusion mechanisms, this work is mainly concerned by the surface, volume and grain boundary diffusions, which are considered to be the most important diffusion routes. Those diffusion mechanisms can be model by using the first Fickss law. This law establish a relation between the matter flow j and the gradient of the chemical potential: where m is the mobility associated with the diffusion path considered. The chemical potential characterize the energy carried by diffusion of a chemical species (here, atoms or vacancies), and its expression depends on region of the particles that is being considered, e.g. surface, volume, grain-boundary, etc. Hence, the matter flux associated with the surface diffusion is expressed by where δ s is the thickness of the layer through which the material diffuses, D s is the surface diffusion coefficient, k the Boltzmann's constant, T > 0 the absolute temperature, γ the surface tension parameter (assumed to be constant here), and finally κ denotes the mean curvature. Note that ∇ s is the surface gradient, and consequently j s is tangential to the grain free surface. Analogously (with obvious notations), the flux associated with the grain boundary diffusion is defined by with σ nn the normal stress. The flux associated with the lattice diffusion is given by where p is the hydrostatic pressure. The mass flux leads to a deposition or removal of mass at the surface and at the grain boundary. This deposition or removal induces a normal velocity which satisfies the following mass conservation equation: v n = −Ω ∇ s ·j for surface and grain boundary routes (5) v n = Ωj · n for lattice route (6) where Ω is the molar volume, and n the unit vector normal to the surface. Previous equations are integrated into a continuum finite element Eulerian approach. The numerical strategy is first presented. Then, the surface and volume diffusion issues are presented. Finally, the bases toward the sintering simulation by grain-boundary diffusion are stated. To conclude, numerical results are shown for two grains and for a granular packing.

Numerical strategy
Initially, the numerical strategy for the simulation of sintering by surface and volume diffusion is presented. The strategy for the simulation of sintering by grain-boundary diffusion is slightly different and will be developed in Section 5. As a first approach only the surface and volume diffusion are considered. For those diffusion routes, the computational domain is formed by only two phases: solid and fluid. Let Ψ be a bounded region, Ψ ⊂ IR d , where d is the spatial dimension. The computational domain Ψ contains two immiscible phases: a set of elastic solids, denoted Ψ s , embedded into a surrounding fluid (air) Ψ f . Furthermore, Ψ s is assumed to be completely contained in Ψ, i.e. Ψ s ∩∂Ψ = ∅. Finally, the interface between the two phases is denoted by Γ sf = ∂Ψ s ∩ ∂Ψ f (see Figure 1).

Fig. 1: Computational domain Ψ
A Level-Set method is used to capture the interface Γ sf . This Eulerian approach consists in embedding the surface in a higher dimensional function, the so called the level-set function. The level-set approach used in this work was proposed by Coupez et al. in [2]. The basic idea of this method is to solve both, the level-set transport and the reinitialization equations, in a single step. The level-set function is initialized by the following expression: for any point (node) x ∈ Ψ, where dist(x, Γ sf ) is the distance from the point x to the interface. Consequently, the interface is described by the zero isosurface of ϕ: Assuming the level-set function ϕ known at a time t, the velocity v n can then be computed by Equation (5), and the level-set function is subsequently transported from t to the time t + ∆t by soling the advection equation: with v = v n n. This equation is solved by a finite element method: Ψ is discretized by a simplex mesh, while ϕ is approximated by a continuous and piecewise linear function. The formulation is then stabilized by a SUPG technique.

Surface diffusion
It is important to note the following metric properties of the level-set function. The values of the levelset function provide the distance from a point x to the interface, the first-order derivative of ϕ provides the normal to the interface : while, the second-order derivatives allow the mean curvature to be computed: However, since ϕ is piecewise linear, its second derivatives vanish, and Equation (10) have to be considered in a weak sense. Furthermore, Equations (2) and (5) show that the surface diffusion velocity is proportional to the fourth-order derivatives of the level-set function! To overcome this difficulty, a mixed stabilized formulation in curvature κ / surface Laplacian of curvature ∆ s κ has been successfully proposed in [3]. The starting point of this formulation is to set: and v n ∝ ∆ s κ. This first-order development introduces a stabilization term ∆ s κ.

Volume diffusion
Mechanical problem. The numerical approach developed to solve mechanical problem is presented in [4]. The inertia terms and the volume forces are neglected. The momentum conservation can be expressed by the set of Equations:

716
The Current State-of-the-Art on Material Forming where σ denotes the Cauchy stress tensor, p ext an external pressure, and [·] sf the operator of jump through Γ sf . Equation (14) is the Laplace's law which expresses the discontinuity of the stress vector through the interface. The solid particles are assumed to behave as isotropic linear elastic bodies. The Cauchy stress tensor as then the form: in Ψ s , where C is the well-known fourth-order elasticity tensor, and ε = (∇ u + ∇ u T )/2 is the strain tensor (small deformations assumption), and u the displacement. As the mechanical problem is formulated in a mixed velocity -pressure form, we first define the mechanical velocity in the solid part as the time derivative of the displacement: and second, we complete the momentum balance equation (12) by defining the hydrostatic pressure through: with tr(·) the trace operator, and K the bulk modulus. The surrounding medium (the air) is assumed to behave as an incompressible Newtonian low viscous fluid. In fact we do not want to deal with the precise dynamics of this medium: in our model, the role of this fluid is only to transmit the pressure p ext , enforced over the domain boundary, to the grain free surface. Thus, Stokes equations are assumed to describe the fluid motion: where η is the fluid viscosity,ε = (∇ v + ∇ v T )/2, and divergence-free relation expresses the mass conservation when dealing with incompressible fluid. The numerical method used to solve the mechanical problem (12) -(18), is presented in [4]. This problem is formulated in terms of velocity and pressure by summing up the solid and fluid variational formulations and accounting for the Laplace's at the interface. Velocity and pressure fields are both approximated by continuous and piecewise linear fields defined on a simplex mesh. In order to deal with discontinuities at the interface (due to Laplace's law but also to the fluid -solid interaction), and to circumvent the Brezzi -Babuska stability condition, a variational multiscale technique is used to stabilize the finite element formulation. However, an important point has to be highlight. The numerical treatment of the surface tension term is not straightforward within an Eulerian approach as the interface Γ sf is not described explicitly by a set of element boundaries. The Surface Local Reconstruction (SLR) method is used. As the interface Γ sf is given by the zero isosurface of the level-set function, it is possible to locally reconstruct a piecewise linear interface. Once the intersection between an element and the interface Γ sf has been found, the surface tension term can then be explicitly computed over this locally reconstructed surface. All the details concerning the numerical implementation of the  [4]. The surface tension over the interface induces a jump of the normal stress in the normal direction. For this reason, the pressure field also presents a jump at the interface between the particles and the surrounding medium. Figure 2 shows the jump of the pressure field atΓ sf as a result of the surface tension present over this interface.
Flow and velocity computation. Equations (2) and (5) show that the lattice flux diffusion depends on the gradient of the pressure computed by solving previous mechanical problem. However, as demonstrated in Figure 2, since one single field is used to describe both fluid and solid pressure fields, this field exhibit a discontinuity at the interface, and consequently its gradient can not be computed in a straightforward way over this interface. This issue is overcome by computing the pressure gradient only inside the grains, up to a distance λ from the interface (see Figure 3). The normal pressure gradient (or normal velocity) computed over {ϕ = −λ} is subsequently projected over the interface {ϕ = 0} by solving: Fig. 3: The normal velocity v n ∝ ∇ p · n is directly computed from the gradient of the pressure in the region ϕ < −λ (darker grey region) and v n is convected in the outward normal direction in the region ϕ > −λ (lighter grey region).

718
The Current State-of-the-Art on Material Forming

Grain boundary diffusion
According to Equations (3) and (5), the velocity associated to this diffusion path is function of the second derivative of the normal stress. The first step toward the sintering simulation by grain boundary is to develop a strategy able to solve the mechanical problem by taking into account, in addition to the surface tension at the free surface of the particles, the surface tension present at the grain boundary. This grain boundary surface tension condition can be written as follows: [σn] gb = γ gb κn (21) over the grain boundary Γ gb , and has to be added to the mechanical problem (12)-(14). Figure 4 shows the pressure field computed by taking into account the surface tension at the free surface and also the surface tension at the grain boundary. Fig. 4: Pressure field computed by taking into account the surface tension over Γ sf and Γ gb (grain boundary). On the right, the jumps of the pressure across the interfaces Γ sf and Γ gb can be seen.
One of the difficulty of simulating grain boundary diffusion is that we have to consider one level-set function per grain. Another difficulty is that the normal stress computed with our strategy is not smooth enough to obtain a good approximation of its Laplacian. For theses reasons, the numerical strategy allowing to compute the velocity due to grain boundary diffusion is yet to be developed. However, the strategy presented for the mechanical problem with surface tension on the free surface and the grain boundary is an important step toward the sintering simulation by grain boundary diffusion.

Numerical results
All the developments have been implemented in the finite element library CimLib. CimLib is a highly parallel C++ library developed at Center for Material Forming (Mines ParisTech, CNRS UMR 7635) by the team of Professor Coupez [5]. As in [6], the simulations presented in this work have been performed by using a mesh adaptation strategy aiming to obtain accurate results while keeping a``reasonable'' number of mesh elements. The mesh is refined over a narrow band around the surface of the particles, in a anisotropic way by using metric properties of the level-set function. In the following, the growth of the neck between two particles is studied, first considering each diffusion path independently and then by coupling both mechanisms. The strategy developed in this work allows to perform simulation over a particles packing, this will be presented at the end of this section.
Two particles sintering: surface and volume diffusions. The first step toward the simulation of sintering process at the particles scale consists in studying the growth of the neck formed between two particles during the early stages of sintering. There exists different models attempting to predict the evolution of the radius of the neck x as a function of the time. Those models are mainly based on some geometrical hypothesis of the evolution of the neck and the particles. In general those geometrical models can be written as follows: where r is the radius of the particles (supposed to remain constant), A is a constant, t * an dimensionless time, and n a parameter which characterizes the diffusion route. Even if the validity of those models is restricted to neck radius smaller that 30% of the radius of the particles, they represent a very useful tool concerning the validation of the kinetic obtained from more elaborated numerical approaches. Figure 5 presents the evolution of the dimensionless radius x/r as a function of t * during surface diffusion for particles radius between 0.1 and 2.5. The exponent n in Equation (22) ranges between 5 and 7 for surface diffusion, which is in agreement with simulations performed. The dashed line showed in Figure 5 (referred to as``Simulation, 1/7'') corresponds to the fitting curve obtained by least-squares and for comparison purposes, the line corresponding to n = 6 in Equation (22) is also plotted (referred to as``1/6 law''). The same kind of simulation has been performed by volume diffusion. The results are shown in Figure 6, with 4 < n < 5, which is in the range predicted by the geometrical models.
Surface and volume diffusions are coupled to simulate the sintering of two particles by this two diffusion mechanisms. Figure 7 shows the evolution of two particles as they sinter together by surface and volume diffusions. Figure 7a presents the initial state of the particles, they are set to be quasitangent. After 50 time steps, the neck between the particles is approximatively 30% of the particle radius, as it is shown in Figure 7b. The power laws presented can used for values of the neck radius x/r < 0.3, it is interesting to see how the two particles evolve beyond this limit (Figure 7c), the particles can not be considered to remain spherical any longer. The coupling between both diffusion mechanisms is showed in Figure 8, and as expected, the neck grows significantly faster when surface and volume diffusion take place simultaneously.
Diffusion in a granular packing. A sintering simulation by volume diffusion over a more realistic particles packing is presented. A set of 154 particles with radius ranging from 0.0633 µm to 0.0797 µm. The set of particles is embedded into a computational domain given by a cube of side 1.2

720
The Current State-of-the-Art on Material Forming  µm. In order to obtain accurate results, the element size is refined near to the surface of the particles. Figure 9a. shows the initial particle packing as well as a cut of the refined mesh that is made up of approximatively 2 millions nodes and 11 millions tetrahedral elements. The evolution of the structure is showed in Figures 9a to 9d. In the initial geometry (9a) particles are set to be quasi-tangent. As the volume diffusion takes place, the necks between the particles grow up to a point (Figure 9d) where the particles can not be distinguish any more. As neither the mass or the density of the particles change during the sintering process, the volume of the particles must remain constant. Considering the simulation showed in Figure 9, the change of total volume of the particles after 200 time steps is of 0.12%, which is negligible. This simulation involves 200 time steps and has been performed in 100 h by using a parallel computing strategy on 24 cores. Fig. 9: Changes occurring in a granular packing by volume diffusion

Conclusion
A numerical approach for the sintering simulation by surface and volume diffusion has been presented. The level-set Eulerian framework adopted allows to take into account all the topological changes that can appear during the sintering. The level-set method also provides very useful tools for the computation of the surface and volume diffusion and at the same time, it allows to handle 2D and 3D simulations with ease. The mesh adaptation strategy used allows to perform simulations while keeping a``reasonable'' number of nodes and elements and is specially useful when dealing with 3D simulations as the number of mesh element is dramatically increased compared with a 2D simulation. The results concerting the sintering of two particles by surface and volume diffusion has been compared with analytical power laws. The comparison showed that both approaches (surface and volume diffusion) lead to a good simulation of the kinetics of each process. Furthermore, specifically regarding the surface diffusion, the results of the evolution of the neck between two particles is in agreement with the predictions obtained from Equation (22). The simulation of sintering by coupled surface and volume diffusions has been presented. Even if the results of this coupling seem to be in agreement with the kinetics of the process, it is very difficult to validate the results as there is not experimental data. Finally the capabilities of the numerical approach has been demonstrated by performing a simulation of the sintering of more realistic particle packing in 3D.

722
The Current State-of-the-Art on Material Forming A large amount of work is yet to be done, mainly concerning the grain boundary diffusion. But the bases for the treatment of the diffusion path has been presented as the approach developed to solve the mechanical problem can also be used to take into account the surface tension at the grain boundary. The mechanical problem has been solved by considering the particles to respond as a linear elastic isotropic material, it would be very interesting to study how different the structure evolves if another material model is considered or if the mechanical properties of each particle are different. The global goal of this work is to simulate the sintering process of fully 3D granular packing, which would open several perspectives regarding the computational materials design.