# Direct 3D simulation of powder sintering by surface and volume diffusion

Daniel Pino Muñoz, Julien Bruchon, Sylvain Drapier, François Valdivieso

► To cite this version:

Daniel Pino Muñoz, Julien Bruchon, Sylvain Drapier, François Valdivieso. Direct 3D simulation of powder sintering by surface and volume diffusion. Key Engineering Materials, 2013, 554-557, pp.714-723. 10.4028/www.scientific.net/KEM.554-557.714. emse-00962456

HAL Id: emse-00962456

https://hal-emse.ccsd.cnrs.fr/emse-00962456

Submitted on 24 Aug 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Direct 3D simulation of powder sintering by surface and volume diffusion

Daniel PINO MUÑOZ¹,ª, Julien BRUCHON¹,ᵇ, Sylvain DRAPIER¹,ᶜ and François VALDIVIESO¹,ᵈ

LCG UMR CNRS 5146, Centre SMS, Ecole Nationale Supérieure des Mines de Saint-Etienne,
158 cours Fauriel, 42023 Saint-Etienne cedex 2, France

ᵃpino@emse.fr, ᵇbruchon@emse.fr, ᶜdrapier@emse.fr, ᵈvaldivie@emse.fr

**Keywords:** solid-state sintering, surface diffusion, volume diffusion, surface tension, level-set.

**Abstract.** Within the general context of solid-state sintering process, this work presents a numerical modelling approach, at the grain scale, of ceramic grain packing consolidation. Typically, the sintering process triggers several matter diffusion routes that are thermally activated: surface, grain boundary and volume diffusions. Including this physics into a high-performance computing framework would permit to investigate and to track the changes occurring into a granular packing during sintering. In performing this kind of simulations, one will face several challenges: the strong topological changes appear during sintering simulation at the grains scale, the evolution of the structure is mainly driven by the surface tension phenomena through the Laplace's law, and the mechanical properties of the grains could, possibly, be different. The proposed numerical simulations are carried out within an Eulerian Finite Element framework and the Level-Set method is used to cope with changes in the microstructure. The results obtained with this numerical strategy are compared with success to the usual geometrical models.

## Introduction

Sintering has become 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 Fick's law. This law establish a relation between the matter flow j and the gradient of the chemical potential:

$$ \mathbf{j} = -m \nabla \mu \qquad (1) $$

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

$$ \mathbf{j}_s = - \frac{\delta_s D_s}{kT} \nabla_s (\gamma \kappa) \qquad (2) $$

where δs is the thickness of the layer through which the material diffuses, Ds 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 $\mathbf{j}_s$ is tangential to the grain free surface. Analogously (with obvious notations), the flux associated with the grain boundary diffusion is defined by

$$ \mathbf{j}_{gb} = - \frac{\delta_{gb} D_{gb}}{kT} \nabla_s \sigma_{nn} \qquad (3) $$

with $\sigma_{nn}$ the normal stress. The flux associated with the lattice diffusion is given by

$$ \mathbf{j}_v = -\frac{D_v}{kT} \nabla p \qquad (4) $$

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 = -\Omega \nabla_s \cdot \mathbf{j} \quad \text{for surface and grain boundary routes} \qquad (5) $$

$$ v_n = \Omega \mathbf{j} \cdot \mathbf{n} \quad \text{for lattice route} \qquad (6) $$

where $\Omega$ is the molar volume, and $\mathbf{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 $\Psi$ be a bounded region, $\Psi \subset IR^d$, where $d$ is the spatial dimension. The computational domain $\Psi$ contains two immiscible phases: a set of elastic solids, denoted $\Psi_s$, embedded into a surrounding fluid (air) $\Psi_f$. Furthermore, $\Psi_s$ is assumed to be completely contained in $\Psi$, i.e. $\Psi_s \cap \partial\Psi = \emptyset$. Finally, the interface between the two phases is denoted by $\Gamma_{sf} = \partial\Psi_s \cap \partial\Psi_f$ (see Figure 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1: Computational domain $\Psi$</figcaption></figure>

A Level-Set method is used to capture the interface $\Gamma_{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:

$$ \phi(\mathbf{x}, t = 0) = \begin{cases} -\mathrm{dist}(\mathbf{x}, \Gamma_{sf}) & \text{if } \mathbf{x} \in \Psi_s \\ +\mathrm{dist}(\mathbf{x}, \Gamma_{sf}) & \text{if } \mathbf{x} \in \Psi_f \end{cases} \qquad (7) $$

for any point (node) $\mathbf{x} \in \Psi$, where $\mathrm{dist}(\mathbf{x}, \Gamma_{sf})$ is the distance from the point $\mathbf{x}$ to the interface. Consequently, the interface is described by the zero isosurface of $\phi$: $\Gamma_{sf} = \{\mathbf{x} \in \Psi : \phi(\mathbf{x}, t) = 0\}$. Assuming the level-set function $\phi$ 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 + \Delta t$ by solving the advection equation:

$$ \frac{\partial \phi}{\partial t} + \mathbf{v} \cdot \nabla \phi = 0 \qquad (8) $$

with $\mathbf{v} = v_n \mathbf{n}$. This equation is solved by a finite element method: $\Psi$ is discretized by a simplex mesh, while $\phi$ 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 level-set function provide the distance from a point $\mathbf{x}$ to the interface, the first-order derivative of $\phi$ provides the normal to the interface :

$$ \mathbf{n} = \frac{\nabla \phi}{\| \nabla \phi \|} \qquad (9) $$

while, the second-order derivatives allow the mean curvature to be computed:

$$ \kappa = \nabla \cdot \mathbf{n} = \nabla \cdot \frac{\nabla \phi}{\| \nabla \phi \|} \qquad (10) $$

However, since $\phi$ 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 $\kappa$ / surface Laplacian of curvature $\Delta_s \kappa$ has been successfully proposed in [3]. The starting point of this formulation is to set:

$$ \kappa^t = \nabla \cdot \frac{\nabla \phi^{t+1/2}}{\| \nabla \phi^{t+1/2} \|} \qquad (11) $$

with $\phi^{t+1/2} = \phi^t + \Delta t \partial \phi^t / \partial t = \phi^t - \Delta t v_n \| \nabla \phi^t \|$, and $v_n \propto \Delta_s \kappa$. This first-order development introduces a stabilization term $\Delta_s \kappa$.

## 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:
$$ \nabla \cdot \underline{\underline{\sigma}} = 0 \text{ in } \Psi \qquad (12) $$

$$ \underline{\underline{\sigma}}\mathbf{n} = -p_{ext}\mathbf{n} \text{ on } \partial\Psi \qquad (13) $$

$$ [\underline{\underline{\sigma}}\mathbf{n}]_{sf} = \gamma\kappa\mathbf{n} \text{ at } \Gamma_{sf} \qquad (14) $$

where $\underline{\underline{\sigma}}$ denotes the Cauchy stress tensor, $p_{ext}$ an external pressure, and $[\cdot]_{sf}$ the operator of jump through $\Gamma_{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:

$$ \underline{\underline{\sigma}}_s = \underline{\underline{C}} : \underline{\underline{\varepsilon}}} \qquad (15) $$

in $\Psi_s$, where $\underline{\underline{C}}$ is the well-known fourth-order elasticity tensor, and $\underline{\underline{\varepsilon}} = (\nabla \mathbf{u} + \nabla \mathbf{u}^T)/2$ is the strain tensor (small deformations assumption), and $\mathbf{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:

$$ \mathbf{v}^t = \frac{\mathbf{u}^t - \mathbf{u}^{t-\Delta t}}{\Delta t} \qquad (16) $$

and second, we complete the momentum balance equation (12) by defining the hydrostatic pressure through:

$$ p = -K \operatorname{tr} \underline{\underline{\varepsilon}} = -K \nabla \cdot \mathbf{u} \qquad (17) $$

with $\operatorname{tr}(\cdot)$ 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:

$$ \underline{\underline{\sigma}}_f = 2\eta\dot{\underline{\underline{\varepsilon}}} - pI \qquad (18) $$

$$ \nabla \cdot \mathbf{v} = 0 \qquad (19) $$

where $\eta$ is the fluid viscosity, $\dot{\underline{\underline{\varepsilon}}} = (\nabla \mathbf{v} + \nabla \mathbf{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 highlighted. The numerical treatment of the surface tension term is not straightforward within an Eulerian approach as the interface $\Gamma_{sf}$ is not described explicitly by a set of element boundaries. The Surface Local Reconstruction (SLR) method is used. As the interface $\Gamma_{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 $\Gamma_{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
<figure><img src="image_3.png" /><figcaption>Fig. 2: Jump of pressure at the fluid - solid interface</figcaption></figure>

method are presented in [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 $\Gamma_{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 $\lambda$ from the interface (see Figure 3). The normal pressure gradient (or normal velocity) computed over $\{\phi = -\lambda\}$ is subsequently projected over the interface $\{\phi = 0\}$ by solving:

$$v_n \propto \nabla p \cdot \mathbf{n} \quad \text{if} \quad \phi < -\lambda \\ \nabla v_n \cdot \mathbf{n} = 0 \quad \text{if} \quad \phi > -\lambda \qquad (20)$$

<figure><img src="image_4.png" /><figcaption>Fig. 3: The normal velocity $v_n \propto \nabla p \cdot \mathbf{n}$ is directly computed from the gradient of the pressure in the region $\phi < -\lambda$ (darker grey region) and $v_n$ is convected in the outward normal direction in the region $\phi > -\lambda$ (lighter grey region).</figcaption></figure>
## 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:

$$[\underline{\underline{\sigma}}\mathbf{n}]_{gb} = \gamma_{gb}\kappa\mathbf{n} \qquad (21)$$

over the grain boundary $\Gamma_{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.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Pressure field computed by taking into account the surface tension over $\Gamma_{sf}$ and $\Gamma_{gb}$ (grain boundary). On the right, the jumps of the pressure across the interfaces $\Gamma_{sf}$ and $\Gamma_{gb}$ can be seen.</figcaption></figure>

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 these 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:

$$ \left( \frac{x(t)}{r} \right)^n = At^* \qquad (22) $$

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'').

<figure><img src="image_6.png" /><figcaption>Fig. 5: Growth by surface diffusion of the dimensionless neck radius $x/r$ over the dimensionless time $t^*$ for different values of $r$ [3]</figcaption></figure>

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 quasi-tangent. 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
<figure><img src="image_7.png" /><figcaption>Fig. 6: Growth by volume diffusion of the dimensionless neck radius x/r over the dimensionless time $t^*$ for different values of r [6]</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Two particles sintering by coupled surface and volume diffusion at: a. $t = 0$, b. $t = 50\Delta t$ and c. $t = 100\Delta t$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8: Geometry of the neck between the particles after 65 time steps</figcaption></figure>

µ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 distinguished 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.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Changes occurring in a granular packing by volume diffusion</figcaption></figure>

## 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 concerning 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.
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.

References

[1] C. Herring, The physics of powder metallurgy, Vol. 27(2) (1951), p.143-179.

[2] L. Ville, L. Silva, T. Coupez, Int. J. Numer. Meth. Fluids, Vol. 66(3) (2011), p.324-344.

[3] J. Bruchon, S. Drapier, F. Valdivieso, International Journal for Numerical Methods in Engineering, Vol. 86(7) (2011), p.845-861.

[4] D. Pino Muñoz, J. Bruchon, S. Drapier, F. Valdivieso, International Journal for Numerical Methods in Engineering, (2012), in press.

[5] H. Digonnet, L. Silva, T. Coupez, in: *AIP Conference Proceedings*, Vol. 908 (2007), p.269.

[6] J. Bruchon, D. Pino Muñoz, F. Valdivieso, S. Drapier, Journal of the American Ceramic Society, Vol. 95 (8) (2012), p. 2398-2405.