# Finite Element simulation for Si3N4 heating and sintering

Thomas Grippi, Stéphanie Béhar-Lafenetre, Holger Friedrich, Daniel Haas, Uwe Schenderlein, Sylvain Marinel, Charles Manière

► To cite this version:

Thomas Grippi, Stéphanie Béhar-Lafenetre, Holger Friedrich, Daniel Haas, Uwe Schenderlein, et al. Finite Element simulation for Si3N4 heating and sintering. Journal of the European Ceramic Society, 2024, pp.116703. 10.1016/j.jeurceramsoc.2024.116703. hal-04634611

HAL Id: hal-04634611

https://normandie-univ.hal.science/hal-04634611v1

Submitted on 4 Jul 2024

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
# Finite Element simulation for Si₃N₄ heating and sintering

Thomas Grippi¹,²*, Stéphanie Béhar-Lafenetre¹, Holger Friedrich³, Daniel Haas⁴, Uwe Schenderlein⁴, Sylvain Marinel², Charles Manière²

1. Thales Alenia Space, 06150, Cannes, France

2. CRISMAT, Normandie Univ, ENSICAEN, UNICAEN, CNRS, 14000, Caen, France

3. Fraunhofer ISC, Center for High Temperature Materials and Design HTL, Bayreuth

4. QSIL Ingenieurkeramic Rauenstein, Germany

**1. Keywords:** Silicon nitride, Si₃N₄, Gas pressure sintering, Densification, Sintering modeling, Finite Element Method

**2. Abstract**

The fabrication of large silicon nitride ceramic components is intricate and demands expertise in gas pressure liquid phase sintering (GPS-LPS). The forefront technology of finite element (FEM) thermo-mechanical modeling plays a key role in sizing aerospace components that necessitate high precision and mechanical strength. A comprehensive sintering simulation should encompass not only densification and grain growth models but also consider the heating environment, accounting for gas convection, conduction, and surface-to-surface radiation. Moreover, silicon nitride sintering simulation addresses the intricate behavior, including the swelling phenomenon during the transitions from intermediate to final stage sintering. The challenge of significant slumping in hollow parts, resulting in costly rejection, is a primary concern for industrial companies when dealing with low-viscosity liquid phase sintering. Validation of the model is achieved through the computation of simple cylindrical cases, while a meticulous comparison between simulations and sintering experiments on complex hollow bodies and deflection bars facilitates precise adjustment of the shear viscosity theoretical parameter derived from the Skorohod-Olevsky continuum theory of sintering. The findings from these simple cases guide the simulation of a complex geometry representative of typical aerospace components, determining the ideal sintering configuration to mitigate slumping issues.

**Graphical abstract**

Finite Element simulation for Si₃N₄ heating and sintering

<figure><img src="image_2.png" /></figure>
### 3. Nomenclature

$\theta$ Porosity

$\dot{\theta}$ Porosity elimination rate (s⁻¹)

$\underline{\sigma}$ Stress tensor (N.m⁻²)

$\dot{\underline{\epsilon}}$ Strain rate tensor (s⁻¹)

$\dot{\epsilon}$ Trace of the strain rate tensor (s⁻¹)

$\varphi$ Shear modulus

$\psi$ Bulk modulus

$P_I$ Sintering stress (Pa)

$\mathbb{I}$ Identity tensor

$\alpha$ Surface energy (J.m⁻²)

$r$ Mean particle radius (m)

$\eta$ Material viscosity (Pa.s)

$\eta_0$ Viscosity pre-exponential factor (Pa.s)

$\gamma$ Bulk modulus exponent

$G$ grain mean diameter (m)

$Q$ Viscosity activation energy (J.mol⁻¹)

$Q_G$ Grain growth model activation energy (J.mol⁻¹)

$k_0$ Grain growth model pre-exponential factor

$R$ Gas constant 8.314 (J.mol⁻¹.K⁻¹)

$\Phi_c$ Convective flux

$\Phi_r$ Radiative flux

$Q$ Viscosity activation energy (J.mol⁻¹)

$R$ Gas constant 8.314 (J.mol⁻¹.K⁻¹)

$T$ Temperature (K)

$T_{ambiant}$ Room temperature (K)

$Q_{PID}$ Regulated heat source (J)

$K_p', K_I', K_D'$ proportional, integrated derivate coefficient

$G_p$ effective shear viscosity (Pa.s)

$K_p$ effective bulk viscosity (Pa.s)

### 4. Introduction

Gas pressure sintering is the most convenient method for the sintering of large complex parts made of silicon nitride.[1,2] This conventional process requires a thermal treatment performed in a gas pressure furnace. In such furnace, the heat is provided by high temperature resistive heating elements, usually graphite [3]. By the means of thermal conduction, convection and radiation, the heat is transmitted to the sample through the gaseous atmosphere[4,5]. The gas pressure is ensured by the $N_2$ to avoid the decomposition of $Si_3N_4$ at high temperature (~1800°C)[6–8]. To reach fully dense parts and due to the covalent nature of Si-N bounding, the silicon nitride ceramics are usually sintered through a liquid phase sintering process, with sintering additives ($Y_2O_3$, MgO, $Al_2O_3$) [9–11]. In the case of large and complex parts, as for
the manufacturing of aerospace parts, this process can imply high distortions due to the very low equivalent viscosities at high temperature resulting from the porous specimen with liquid phase [12,13]. A simulation of process coupling the heating environment and the densification aspects should be a cutting-edge tool for the manufacturing of such parts as it would replace the long and expensive trial and error tests by virtual sintering tests.

The best simulation tools may be capable for predicting heat transfer by convection, surface to surface radiation, conduction in solid parts and for the sintering phenomena, the specimen densification and microstructure development. To make the simulation as close as possible to the experiment even the proportional integral derivative thermal regulation of the furnace should be considered. Numerous industrial furnace simulations can be found in the literature. These simulations are based either on finite element or on finite volume methods [14–21]. For this study our previous study ref[14] is used as a basis as it used a fully coupled fluid dynamics thermo-mechanical model to predict the uniformity of heating of a large gear shape specimen. The simulation tool should also relate the temperatures and their potential heterogeneities with the sintering. Continuum theory of sintering should be used for this endeavor as this kind of model can be integrated in a finite element code. Numerous models exist [22] from Skorohod-Olevsky [23], Abouaf [24,25], Kraft and Riedel [26] to viscous poisons ratio form [27,28]. In this context, the sintering of silicon nitride has a particularly complex behavior implying a high grain growth sensitivity [29] and a complex swelling phenomenon that have been previously study and modeled analytically [30]. The latter are used is the present study to simulate complex shapes sintering.

The objective of the latter is to depicts the use and the adaptation of a similar coupled thermo-mechanical model for silicon nitride sintering, adding the problematics of gas pressure liquid phase sintering and gravity implied distortions.

In this article, a verification of the analytical sintering model is presented, the application of the model to complex shape aerospace parts is the ultimate goal of the study. The authors first aim to go further from pure analytical model through simulations and experiments. The sintering of
a simple cylinder is studied for the coupling of a mechanical sintering model (with grain growth and densification model [29,30]) and a heat transfer model. Above mentioned grain growth and densification models are both part of the group's global study on silicon nitride sintering for the production of large and complex shapes for aerospace structures. The latter coupled model is applied to more complex geometries to compute and adjust shear modulus and thus, deflection behavior. Additionally, in preparation for applying the model to large, complex aerospace parts, the authors intend to improve the accuracy of sintering modeling by integrating it with other physical phenomena. For large furnaces, it is crucial to monitor the sintering environment, temperature gradients, and gas recirculation to ensure homogeneous sintering. The ultimate study focuses on the application of this model to a real industrial part. The final objective of the depicted study is the identification of a most suitable industrial sintering configuration with the particular case of an industrial component.

# 5. Experiments and Method

## 5.1. *Equipment*

The gas pressure furnace to be simulated is TOM_Metal from the Fraunhofer-center HTL. This gas pressure furnace is a vertically placed cylinder with radial graphite heaters. The usable volume of the furnace is of 200mm diameter and 400mm height. The reference material studied is a compacted and debinded form of “gas pressured sintered grade” of silicon nitride produced by QSIL Ingenieurkeramik, depicted in the authors’ previous publication[30]. The studied samples are debinded cold isostatic pressed cylinders of ø10×10mm. These are placed into 2 interlocked graphite crucibles one with windows for the camera (TOM) dilatometry devices[31] that follows the specimen sintering and another control sample in a fully closed crucible. The gas used for the pressure sintering, to avoid decomposition of Si₃N₄ at 1800°C, is N₂ [6]. In previous studies[29,30], the samples from this same grade were used for the parameters
extraction and the conception of a sintering model have been compacted and debinded by QSIL. Studied samples for the FEM simulation are ø10x10mm cylinders for verifications of the coupled model. A range of hollow cylinders (ø50x50mm, ø40x40mm and ø30x30mm both with a 4mm thickness), placed horizontally and 40x10mm deflection bars with 2 mm/4 mm thickness were prepared to study gravity implied distortions. A more complex geometry for an industrial case has also been computed in a final stage.

## 5.2. Physics simulation:

As the simulation study requires the coupling of several physics, it has been conducted on the software COMSOL Multiphysics. In this simulation, the sintering model (densification and the grain growth) is coupled to the heat transfer model. In the present case, heat transfer contributions are thermal radiations, implemented from heat equations, and convections implemented from heat equations and Navier-Stokes equations[4] (as detailed in Figure 1).

<figure><img src="image_3.png" /><figcaption>Figure 1 Schematic of the multiphysic modeling strategy and the shared parameters</figcaption></figure>

The densification model is based on Skorohod-Olevsky continuum theory of sintering constitutive equation [23,32,33].:

$$ \underline{\sigma} = \frac{\sigma(w)}{w} \left( \varphi \dot{\underline{\epsilon}} + \left( \psi - \frac{1}{3} \varphi \right) \dot{\underline{e}} \right) + (P_l - P_s) \dot{\underline{i}} \quad (1) $$
Where $\sigma(w)$ is the equivalent stress and $w$ the equivalent strain rate, $\dot{\epsilon}$ is the strain rate tensor, $\ddot{\mathbf{i}}$ is the identity tensor, $\varphi$ and $\psi$ the shear and bulk moduli, the invariant $\dot{e} = \dot{\epsilon}_x + \dot{\epsilon}_y + \dot{\epsilon}_z$, Pl the capillary stress, $P_s$ the swelling pressure opposing to the densification. The mechanical model has been adapted for silicon nitride liquid phase sintering under gas pressure, detailed in a previous article [30]. This previous article details the extraction of sintering parameters for the conception of a model that takes into account both the granular growth and the powder densification.

In our particular case, the slight gas over pressure sintering can be assimilated to free sintering, the $\frac{\sigma(w)}{w}$ term thus becomes $2\eta$ for linear viscous behavior:

$$ \underline{\sigma} = 2\eta \left( \varphi \underline{\dot{\epsilon}} + \left( \psi - \frac{1}{3}\varphi \right) \dot{\epsilon} \ddot{\mathbf{i}} \right) + P_s \ddot{\mathbf{i}} - P_s \ddot{\mathbf{i}} \quad (2) $$

With the Skorohod theoretical shear and bulk moduli expressed as:

$$ \varphi = (1 - \theta)^2 \quad (3) $$

$$ \psi = \frac{2(1-\theta)^3}{3\theta} \quad (4) $$

The expression of the bulk moduli for the present case has been slightly modified based on experimental measurements and derivation methods detailed in a previous study [30].

This more comprehensive expression takes into consideration the porosity at the initial stage of sintering by the mean of this critical porosity parameter $\theta_c$ [24,34]. The initial critical porosity term accounts for the porosity at the initial stage of sintering, providing a more realistic expression than the theoretical Skorohod' expression. In the latter expression, this term is set to unity, implying 100% porosity at initial stage of sintering which can lead to modeling bias.

$$ \psi = \frac{2(\theta_c - \theta)^\gamma}{3\theta} \quad (5) $$

Based on the linear diffusional creep mechanisms, the viscosity grain size sensitivity and temperature dependence can be defined by the expression [35]:

$$ 2\eta(G, T) = \left(\frac{G}{G_0}\right)^m \eta_0 T \exp\left(\frac{Q}{RT}\right) \quad (6) $$
In the present case of viscous flow induced by liquid phase sintering, the term m (viscosity grain size exponent) is set to 0[36]. The last expression is then reduced to:

$$2\eta(T) = \eta_0 T \exp\left(\frac{Q}{RT}\right) \qquad (7)$$

From this linearization and the adjustment of the bulk modulus the following parameters were obtained [30] Q=340kJ/mol, $\eta_0$ = 0.021, $\theta_c$ =0.6 and $\gamma$=1.24.

The expression of the sintering stress P₁ can be developed as follows, evolving as a function of the density and the grain size.

$$P_l = \frac{3\alpha}{G} (1 - \theta)^2 \qquad (8)$$

The evolution of the grain size is also considered with the model depicted in the article[29] and is expressed as:

$$\dot{G} = \frac{K(T)}{G^P} \quad \text{with} \quad K(T) = k_0 \exp\left(-\frac{Q_G}{RT}\right) \qquad (9)$$

The latter grain growth model encompasses 2 different kinetic regimes as a function of the heating and the temperature holding. This transition is held with a mathematical switch between the couple of parameters for the heating $Q_{G1}$=262kJ/mol and $k_{01}$=1,009E-15 m$^{1+p}$.s$^{-1}$ and the couple of parameters for the temperature holding $k_{02}$ = 7,41575E-13 m$^{1+p}$.s$^{-1}$ and $Q_{G2}$ = 399kJ/mol:

$$\left\{ \begin{array}{l} f_0 = \frac{1}{1+e^{-(t-t_0)}} \\ Q_G = Q_{G1}(1-f_0) + Q_{G2}(f_0) \\ k_0 = k_{01}(1-f_0) + k_{02}(f_0) \end{array} \right. \qquad (10)$$

As developed in the latter article[30], observations and optical dilatometries on the studied silicon nitride based ceramics displayed important swelling that has been associated to 2 possible mechanisms. Either non-equiaxial grain growth (linked to the $\alpha/\beta$-Si$_3$N$_4$ allotropic transition)[6,37] and pore coalescence [36] have been explored generating a swelling phenomenon. Indeed, the present model takes into account the gas trapped in the pores and the resulting pressure that could slow the densification. In this case the general model is taking into account more parameters. With Equation (3), the equivalent analytic model is the following
$$ \dot{\theta} = \frac{(-Pl+Ps+\sigma_r)(1-\theta)}{\eta_0 T \exp\left(\frac{Q}{RT}\right) \psi(\theta)} \quad (11) $$

In the latter analytical model, the swelling pressure ($P_s$) is considered, representing the pressure generated from swelling that counters both the sintering pressure ($P_t$) and the externally applied gas pressure ($\sigma_r$) which aids in eliminating porosity in the final stage (approximately 5 bar). The latter expression includes the Equations 2,3,5,7,8,10. The swelling pressure which is the last remaining parameter to identify has been assessed using the following parameters: the grain size evolution (grain growth model), the identified viscosity and moduli and the experimentally imposed gas pressure evolution ($\sigma_r$). Using the data on shrinkage during the final stage of sintering, we can determine pore swelling by rearranging the model Equation (11) to isolate the $Ps/\alpha$ ratio for calculation.:

$$ \frac{Ps}{\alpha} = \frac{\dot{\theta} \frac{\eta_0 T \exp\left(\frac{Q}{RT}\right) \psi(\theta)}{(1-\theta)}}{6(1-\theta)^2} - \frac{\sigma_r}{\alpha} \quad (12) $$

The whole furnace of the Fraunhofer-center HTL, TOM_Metal, has been reproduced in a 2D-axisymmetric model (Figure 2). In the enclosure of the furnace, crucibles, graphite heaters and studied samples in the 2 configurations have been represented [5,14]. The witness sample configuration (Figure 2 a and b) experimentally aided in comparing the sample in an open crucible to one in closed crucible conditions. For the purpose of optical dilatometry, the main specimen was placed in an open crucible, which could be more subject to environmental and gaseous perturbations than parts industrially sintered in a closed crucible.

Insulating parts of the furnace have been also computed: graphite fiber insulator and water cooling on top and bottom covers: The convective and radiative cooling fluxes on the external surfaces is expressed by:

$$ \Phi_c = h_{conv}(T_{ambiant} - T); \quad \Phi_r = \sigma_s \varepsilon (T^4 - T_{ambiant}^4) \quad (13) $$

With $h_{conv}$ the convective heat transfer coefficient, $\sigma_s$ the Stefan-Boltzmann's constant, $\varepsilon$ the emissivity, $T_{ambiant}$ the room temperature 25°C and T the temperature on the surface.
The heat source $Q_{PID}$ is regulated with a PID whose expression is [38]:

$$Q_{PID} = K_{P}e(t) + K_{I} \int_{0}^{t} e(t)dt + K_{D} \frac{de(t)}{dt} \quad (14)$$

With coefficient $K_P=400$, $K_I = 20$ and $K_D = 1000$, providing a good regulation.

And In the model, the gas-pressure cycle has been regulated with an outlet, in a similar way as it is done on TOM_Metal. Crucibles, supports, and heating elements have been implemented as graphite parts[3]. As mentioned above, insulators are in graphite fibers and thermal properties have been found in a furnace supplier database [39]. Concerning the material properties of the furnace's frame, a common steel has been chosen to fit to the actual TOM_Metal[40]. Top and bottom covers are also in the same steel and water-cooled (see Figure 2). The gaseous atmosphere of sintering is $N_2$ and silicon nitride samples have been implemented with such properties as thermal conductivity and specific heat capacity given below [23,41–47]:

Thermal conductivity as a function of temperature:

<table><caption>Table 1 Silicon Nitride thermal conductivity</caption><thead><tr><th>T(K)</th><th>Thermal conductivity(W/m.K)</th></tr></thead><tbody><tr><td>301</td><td>79</td></tr><tr><td>551</td><td>52</td></tr><tr><td>840</td><td>39</td></tr><tr><td>1114</td><td>33</td></tr></tbody></table>

Specific heat capacity as a function of temperature

<table><caption>Table 2 Silicon nitride heat capacity</caption><thead><tr><th>T(K)</th><th>Heat capacity(J/kg.K)</th></tr></thead><tbody><tr><td>298</td><td>690</td></tr><tr><td>423</td><td>887</td></tr><tr><td>573</td><td>1030</td></tr><tr><td>873</td><td>1170</td></tr><tr><td>1273</td><td>1270</td></tr></tbody></table>

Other properties that are not considered as temperature dependent in our model, have also been computed in the silicon nitride properties part of the software:
<table><caption>Table 3 Material properties for silicon nitride</caption><thead><tr><td>Property</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>Density</td><td>3100</td><td>kg/m³</td></tr><tr><td>Young's modulus</td><td>325</td><td>GPa</td></tr><tr><td>Poisson's ratio</td><td>0.26</td><td>1</td></tr><tr><td>Coefficient of thermal expansion</td><td>2.3x10<sup>-6</sup></td><td>1/K</td></tr><tr><td>Relative permittivity</td><td>9.7</td><td>1</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>Figure 2 Furnace sintering configurations (a.), geometry of the 2D-Axisymmetrical coupled model (b.) and meshing of the geometry (c.).</figcaption></figure>

Simple ø10x10mm cylinders have been studied for the model verification. Other geometries such as hollow cylinders and deflection bars were sintered to adjust shear modulus from Skorohod model. A comparison between the final experimental deflection values and the simulated geometry with the coupled model. The adjustment of the model to simulate the right deflection is done through variations of the effective shear viscosity ($G_p = \eta(T)\phi$).

The finality of the whole study, is to be able to compute real complex part, in an industrial case.
The first demonstrator to be computed has been provided by Thales Alenia Space (Figure 3).

This simulation will help for the distinction of a most suitable sintering configuration and/or an
acceptable sizing.

<figure><img src="image_5.png" /><figcaption>Figure 3 CAD of the industrial part and its contact area with the support, colored in green for two distinct sintering configuration “cantilever” above and “bridge” below</figcaption></figure>
# 6. Results

## 6.1. Model validation

This sub part verifies the accuracy of the above developed coupled thermo-mechanical model of sintering, with a comparison between experimental densification curves (extracted from optical dilatometries) and simulated densification. As can be observed on Figure 4, the whole furnace enclosure is modeled through a heat-transfer model. With this configuration that differs from the purely mechanical simple cylinder case, the response to the sintering cycle appears to be different from the analytical model detailed in a previous article[30]. Here, the temperature seen by the samples differs from the temperature measured by the pyrometer (Figure 4 a.). The shading crucibles and the distance of the sample from heaters and the pyrometer target are reasonable explanations for this shift. At low temperature, this shift has no effect on the sintering, as it starts to be effective at 1350°C. However, the remaining temperature gap of 10°C to 15°C, between the measured temperature and the temperature seen by the samples, needs to be recomputed for the calculation of the parameter $P_s$, which is very sensible to the this temperature variation. From this swelling pressure adjustment with the simulated temperature, the behavior of the sample is well computed. As can be observed below on Figure 4 d., the simulated behavior approaches in a better way the experimental values. The granular growth model developed in this publication [29] and implemented in this publication[30], can also be computed in the FEM software. This first FEM simulation succeeds in representing the heating and sintering data and fitting the experiment. Adjustments on $P_s$ based on a simulated temperature are required to fit as much as possible to the reality. This is a right basis for the study of the sintering of more complex geometries. The computed thermo-mechanical model shows comprehensive response on these simple cases of 10x10 cylinders. As confirmed experimentally, the latter samples are not subject to significant perturbation such as gravity induced distortions [48], thus, this model basis will be developed in the following subsection to comprehend accurately these phenomena.
<figure><img src="image_6.png" /></figure>

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

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

Figure 4 Temperature cycle, shift between the temperature PID and the temperature in the crucibles (a.), complete modeling of the furnace with heat-furnace and CFD (b.), densification model (c.) and simulated relative density as a function of time (d.).

## 6.2. Shear modulus adjustment

The aim of this section is the study of distortions in more complex parts and adjust parameters of the model to fit to the reality. The case of cantilever bending bars is the first to be modeled, to adjust parameters influencing the deflection during the sintering due to gravity effect, then a case of a hollow cylinder is computed for verification. Bars are placed in a 30 mm cantilever position, as shown in Figure 5. This study has been done in 2 similar steps, 2 bars are placed on a graphite support. One of the bars is weighted with a ø10x10mm debinded cylinder (Figure 5. left) and the other is in a free bending configuration (Figure 5. right). As the 2 mm bars were too thin i.e. the experiment has failed and the mass bar touched the ground, the only considered case for simulation is the 4 mm thickness bars.
<figure><img src="image_9.png" /><figcaption>Figure 5 Bending sintering configuration with 2 mm thickness bar (picture a. at debinded state, picture b. sintered) and 4 mm thickness bar (picture c. at debinded state, picture d. sintered)</figcaption></figure>

The mesh configuration stands out as the first critical parameter for investigation. Figure 6 illustrates the high sensitivity of bar deflection to mesh quality. The linear discretization (blue) is particularly sensitive, showing significantly lower deflection as the mesh transitions from coarse to fine. Conversely, the quadratic discretized mesh (orange) consistently produces similar results regardless of mesh refinement. It's noteworthy that the number of degrees of freedom is directly correlated with computational time. As depicted in Figure 6, selecting a linear mesh necessitates a substantial level of refinement, resulting in extended computational time to approach the stable results obtained with a quadratic mesh. Therefore, simulating deflection or complex sintering distortions using a quadratic mesh is imperative to avoid underestimating the magnitude of distortion.
<figure><img src="image_10.png" /><figcaption>Figure 6 Influence of mesh configuration on the bar deflection test.</figcaption></figure>

A simulation of the bending configuration is performed by applying the developed and checked sintering model to a similar geometry as the experiment (see Figure 7 b.). From this simulation a deflection value can be extracted and compared to the experiment (see Figure 7 a.).

In the analytical model developed in the article [30], based on simple geometry, the shear behavior is intentionally limited to simplify the study and identify the main constitutive parameters. In this original analytical case, the expression of Skorohod's shear modulus is kept in its theoretical form Equation (3). However, while the final shape and density can be predicted with the FEM for simple geometries such as small bulk cylinders (10x10mm), we observed a discrepancy for hollow parts or cantilever bodies where the shear behavior is no longer negligible. Therefore, by making a slight correction to the shear modulus, it is possible to adjust the shear behavior of the part while maintaining its correct densification behavior. This action consists in increasing the high temperature deformation sensitivity without touching the sintering densification behavior identified. The comparison is made to adjust the model and simulate the right deflection is done through variations of the effective shear viscosity ($G_p$= $\eta(T)\varphi$). A factor $c$ is added to the expression of Skorohod's shear modulus to affect the shear viscosity behavior:

$$ \varphi = c * (1 - \theta)^2 \qquad (15) $$

With $\theta$ being the porosity and $c$ a coefficient that equals 1 for Skorohod theoretical law[49].

This approach allows to adjust the distortion sensitivity independently from the previous identified densification behavior ($K_p=2 \eta(T)\psi$). As can be observed on the Figure 8 the common
expression of the shear modulus with $c=1$ does not permit the model to approach the experimental values. Several cases are studied from $c=1$ to $c=2$ and the simulation fits the best the experimental values in both configuration with $c=1.85$.

<figure><img src="image_11.png" /><figcaption>Figure 7 Comparison between post-sintering measured deflection with a 3D scan (a.) and simulated deflection (b.)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 8 Shear modulus $\phi$ adjustment from free deflection and loaded deflection experimental data.</figcaption></figure>
### 6.3. Hollow cylinders

From the adjusted sintering model, a second comparison between experimental hollow
cylinders samples and simulated ones has been performed (Figure 9). The principal objective
of this secondary study is the verification of the adjusted shear modulus on a significantly
different geometry also subjected to gravity. Hollow cylindrical shapes can frequently be
encountered during complex industrial parts design and sizing. This secondary computation
gathers supplementary material for the verification of the model's robustness. The model used
in this study is a 3D model with the same mesh size and refinement as for the bending study.
By keeping the same parameters, and indeed the same φ value as adjusted in the previous
section, the model is able to approach the distortions implied by gravity with an error of 1% for
the Ø50x50. Results between the comparison of sintered samples and simulation of the same
samples for the 3 cases Ø50x50x4mm, Ø40x40x4mm and Ø30x30x4mm can be seen in the Table
4. The model is capable of reproducing the effect of the gravity for all cases. Due to its more
important thickness/diameter ratio, the Ø30x30 cylinder is the less subject to distortions. These
thickness aspects and gravity-implied distortions should be carefully studied during the design
of a large complex hollow part otherwise it could not respect final industrial geometrical
specifications and will be rejected.
<figure><img src="image_13.png" /><figcaption>Figure 9 Experimental sintering case of hollow cylinders (pictures a. and b.), final shape from sintering simulation of 3 cylinders : φ50x50x4mm (picture c. and f.), φ40x40x4mm (picture d. and g.), and φ30x30x4mm (picture e. and h.) displacement field on z axis (Above line), total displacement field(Below line).</figcaption></figure>
<table><caption>Table 4 Gap/error calculation between sintered samples and simulation.</caption><thead><tr><th></th><th>debinded<br>diameter 1 &<br>2 (mm)</th><th>sintered<br>diameter 1<br>(mm)</th><th>sintered<br>diameter 2<br>(mm)</th><th>Gap<br>simulation/reality<br>diameter 1 (%)</th><th>Gap<br>simulation/reality<br>diameter 2 (%)</th></tr></thead><tbody><tr><td>EXP Cyl. 50x50</td><td>50,1</td><td>40,16</td><td>41,18</td><td></td><td></td></tr><tr><td>SIM Cyl. 50x50 with φ=1.85(1-θ)<sup>2</sup></td><td>50,1</td><td>39,72</td><td>41,37</td><td>1,10%</td><td>-0,46%</td></tr></tbody></table><table><thead><tr><th></th><th>debinded<br>diameter 1 &<br>2 (mm)</th><th>sintered<br>diameter 1<br>(mm)</th><th>sintered<br>diameter 2<br>(mm)</th><th>Gap<br>simulation/reality<br>diameter 1 (%)</th><th>Gap<br>simulation/reality<br>diameter 2 (%)</th></tr></thead><tbody><tr><td>EXP Cyl. 40x40</td><td>40,05</td><td>32,37</td><td>32,74</td><td></td><td></td></tr><tr><td>SIM Cyl. 40x40 with φ=1.85(1-θ)<sup>2</sup></td><td>40,05</td><td>32,2</td><td>32,78</td><td>0,53%</td><td>-0,12%</td></tr></tbody></table><table><thead><tr><th></th><th>debinded<br>diameter 1 &<br>2 (mm)</th><th>sintered<br>diameter 1<br>(mm)</th><th>sintered<br>diameter 2<br>(mm)</th><th>Gap<br>simulation/reality<br>diameter 1 (%)</th><th>Gap<br>simulation/reality<br>diameter 2 (%)</th></tr></thead><tbody><tr><td>EXP Cyl. 30x30</td><td>30,02</td><td>24,32</td><td>24,42</td><td></td><td></td></tr><tr><td>SIM Cyl. 30x30 with φ=1.85(1-θ)<sup>2</sup></td><td>30,02</td><td>24,24</td><td>24,4</td><td>0,33%</td><td>0,08%</td></tr></tbody></table>

## 6.4. Complex shape

The last case to be studied is the demonstrator part provided by the industrial partner. A CAD model is used in this case with the same mesh size and refinement applied as for the previous studies. The mechanical model developed is thus applied to the FEM mesh. Two distinct sintering configurations can be studied easily, a “cantilever” configuration and a “bridge” configuration. These differs one from each other by their sintering configuration orientation. In the first case, (“cantilever”) the part is sintered placed horizontally with contacts with the support as shown on Figure 3 above. In the second case, named for the study “bridge”, the same part is rotated by 180° on the y-axis and the contacts with the support are highlighted on the picture Figure 3.

The simulation results, depicted in both "side" and "3D" views, are presented in Figure 10a for the "cantilever" configuration and in Figure 10b for the "bridge" configuration. It is evident that the "cantilever," lacking adequate support at the edge of the curved part, experiences preferential slumping in this region. This can be clearly observed in the "side" views. Figure 10c illustrates the displacement probe, clearly indicating this distortion with a black curve for the "bridge" configuration, which closely resembles the sintering curve, and a significantly higher displacement magnitude observed during the dwell for the "cantilever" case.
However, the most effective means of detecting these distortions remains the comparison
between the FEM-generated final "cantilever" shape and the ideal sintered CAD part without
any distortion (achieved through a homothetic transformation of the original geometry). This
CAD dimensional comparison is presented in Figure 10d. The error bars represent the surface-
to-surface gap between these two shapes, revealing distortions of up to 4mm in magnitude for
the unsupported zones. The central zones and the two tubular zones at the base of the shape are
well supported and closely resemble the ideal shapes represented in (dark translucent shape).

The simulation tool enables the quantification of distortions for CAD design adjustments and facilitates the selection of the most optimal sintering configuration. This outcome underscores the relevance of the simulation tool for liquid-phase sintering processes, such as the sintering of silicon nitride, which are highly sensitive to distortions if a part of the shape lacks proper support. The effectiveness of such a comprehensive tool is very promising for optimizing the sizing and CAD design of large, complex shapes in the industry
<figure><img src="image_14.png" /><figcaption>a. Cantilever support configuration</figcaption></figure>

Max 10.3

Min 0.825

<figure><img src="image_15.png" /><figcaption>b. Bridge support configuration</figcaption></figure>

Min -3.73

Max 0

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

Max 10.3

Min 0.825

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

Min -3.73

<figure><img src="image_18.png" /><figcaption>d. Shape distortions error analysis (CAD) for the cantilever configuration</figcaption></figure>

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

Figure 10 Finite Element Method (FEM) simulations were conducted for sintering complex-shaped objects in both (a) cantilever and (b) bridge configurations. To visualize the magnitude of distortion in both configurations, a displacement probe along the z-axis was placed on the object at the location indicated by the red dot. The resulting displacement is then reported in (c) and represents the sum of the sintering shrinkage displacement and the distortion displacement. In the case of the cantilever configuration, the distortion is significant and can be compared to idealized shapes shown in (d). The colored shape represents the simulated shape, with the color bar indicating the error in millimeters compared to the ideal sintered shape (dark translucent)..
7. Conclusion

Based on a mechanical model of sintering that couples a densification and grain growth model, along with a reliable furnace heat transfer model, it was possible to simulate the gas pressure sintering of simple silicon nitride parts. Gravity-induced distortions, stemming from liquid phase sintering, pose significant challenges in industrial sizing processes. The sintering of silicon nitride necessitates numerous experimental iterations to mitigate distortions that arise during the process. Therefore, the remainder of the study focuses on evaluating the model's ability to replicate sintering deformations in various configurations. The investigation was conducted at three levels: bar deflection, hollow cylinders, and a complex shape.

The initial bar deflection analysis highlights the importance of using a quadratic mesh to avoid underestimating sintering distortions. This test reveals that the dilatometry-based model underestimates sintering distortions, and a correction to the shear modulus can be implemented to address this issue. The resulting model provides an excellent fit for the distortions observed in loaded bar and hollow cylinder experiments.

Furthermore, the examination of a complex shape demonstrates how the simulation tool can be applied to analyze sintering configurations. Depending on the shape and orientation within the furnace, sintering zones may receive varying levels of support, potentially leading to severe distortions if not properly managed. Numerical optimization of the sintering configuration is achievable by integrating CAD design, supports, advanced topology optimization tools, and the simulation tool discussed in this study.
## 8. Acknowledgement

We would like to thank Thales Alenia Space for the PhD subvention of the author (CIFRE 2019/1366). A collaboration between four entities has been needed, QSIL-Ingenieurkeramik (industry) has provided the bound specimens, the advanced gas sintering equipment of the Fraunhofer center HTL has been used for the experimental sintering test, the modeling/simulation study is processed at the CRISMAT laboratory, and the assessment of the results for space applications is done by Thales Alenia Space. The help and support of Waldemar Walschewski is gratefully acknowledged.

## 9. Credit authorship contribution statement

**Thomas Grippo:** Conceptualization, Supervision, Experimentation, Modeling, Writing;
**Stephanie Behar-Lafenetre:** Conceptualization, Supervision, Review & Editing;
**Daniel Haas:** Experimentation, Review & Editing;
**Uwe Schenderlein:** Experimentation, Review & Editing;
**Holger Friedrich:** Conceptualization, Experimentation, Supervision, Review & Editing;
**Sylvain Marinel:** Supervision, Conceptualization, Review & Editing.
**Charles Maniere:** Supervision, Conceptualization, Modeling, Writing, Review & Editing

## 10. References

[1] R.M. German, History of Sintering, Sinter. from Empir. Obs. to Sci. Princ. (2014) 13–40. https://doi.org/10.1016/b978-0-12-401682-8.00002-1.

[2] R.M. German, History of sintering: Empirical phase, Powder Metall. 56 (2013) 117–123. https://doi.org/10.1179/1743290112Y.0000000025.

[3] V. V Kleimenov, V. V Alekseev, Carbon heating elements for high-temperature, (1993) 50–51.

[4] A. Minea, A review on analytical techniques for natural convection investigation in a heated closed enclosure: Case study, Therm. Sci. 19 (2015) 1077–1095. https://doi.org/10.2298/TSCI131027021M.
[5] C. Manière, F. Borie, S. Marinel, Impact of convection and radiation on direct/hybrid heating stability of field assisted sintering, J. Manuf. Process. 56 (2020) 147–157.
https://doi.org/10.1016/j.jmapro.2020.04.075.

[6] F.L. Riley, Silicon Nitride and Related Materials, J. Am. Ceram. Soc. 83 (2004) 245–265.
https://doi.org/10.1111/j.1151-2916.2000.tb01182.x.

[7] S.K. Biswas, F.L. Riley, Gas pressure sintering of silicon nitride - current status, Mater. Chem. Phys. 67 (2001) 175–179.
https://doi.org/10.1016/S0254-0584(00)00436-3.

[8] G. Petzow, M. Herrmann, Silicon Nitride Ceramics, in: Eng. Ceram. Curr. Status Futur. Prospect., 2002: pp. 47–167.
https://doi.org/10.1007/3-540-45623-6_2.

[9] V. V Krasil, V. V Sirota, A.S. Ivanov, L.N. Kozlova, O.A. Luk, INVESTIGATION OF THE STRUCTURE OF Si 3 N 4 -BASED CERAMIC WITH Al 2 O 3 AND Y 2 O 3 ADDITIVES, 71 (2014) 15–17.

[10] F. Zhou, J. Pan, K. Chen, Liquid-phase bonding of silicon nitride ceramics using Y2O3-Al2O3-SiO2-TiO2 mixtures, Mater. Lett. 58 (2004) 1383–1386.
https://doi.org/10.1016/j.matlet.2003.09.032.

[11] M. Mitomo, Pressure sintering of Si3N4, J. Mater. Sci. 11 (1976) 1103–1107.
https://doi.org/10.1007/BF02396645.

[12] E.A. Olevsky, R.M. German, Effect of gravity on dimensional change during sintering—
I. Shrinkage anisotropy, Acta Mater. 48 (2000) 1153–1166.
https://doi.org/10.1016/S1359-6454(99)00368-7.

[13] R.M. German, P. Suri, S.J. Park, Review: Liquid phase sintering, J. Mater. Sci. 44 (2009) 1–39.
https://doi.org/10.1007/s10853-008-3008-0.

[14] C. Manière, T. Zahrah, E.A. Olevsky, Fluid dynamics thermo-mechanical simulation of sintering: Uniformity of temperature and density distributions, Appl. Therm. Eng. 123 (2017) 603–613.
https://doi.org/10.1016/j.applthermaleng.2017.05.116.

[15] R. Prieler, S. Pletzer, S. Thusmer, G. Schwabegger, C. Hochenauer, Numerical simulation of a fire resistance test and prediction of the flue gas leakage using CFD/FEM
coupling, J. Struct. Fire Eng. (2023). https://doi.org/10.1108/JSFE-01-2023-0011.

[16] T. Detemmerman, F. Froment, Three dimensional coupled simulation of furnaces and reactor tubes for the thermal cracking of hydrocarbons, Rev. l'Institute Fr. Du Pet. 53 (1998) 181-194. https://doi.org/10.2516/ogst:1998017.

[17] J. Arroyo, L. Pérez, V. Cuervo-Piñera, CFD Modeling and Validation of Blast Furnace Gas/Natural Gas Mixture Combustion in an Experimental Industrial Furnace, Processes. 11 (2023) 1-24. https://doi.org/10.3390/pr11020332.

[18] J.H. Chang, J. Oh, H. Lee, Development of a roller hearth furnace simulation model and performance investigation, Int. J. Heat Mass Transf. 160 (2020) 120222. https://doi.org/10.1016/j.ijheatmasstransfer.2020.120222.

[19] G. Daurer, J. Raič, M. Demuth, C. Gaber, C. Hochenauer, Detailed comparison of physical fining methods in an industrial glass melting furnace using coupled CFD simulations, Appl. Therm. Eng. 232 (2023) 121022. https://doi.org/10.1016/j.applthermaleng.2023.121022.

[20] C. Yang, G. Liu, C. Chen, Y. Hou, M. Xu, Y. Zhang, Numerical simulation of temperature fields in a three-dimensional SiC crystal growth furnace with axisymmetric and spiral coils, Appl. Sci. 8 (2018). https://doi.org/10.3390/app8050705.

[21] E. Hachem, G. Jannoun, J. Veysset, M. Henri, R. Pierrot, I. Poitrault, E. Massoni, T. Coupez, Modeling of heat transfer and turbulent flows inside industrial furnaces, Simul. Model. Pract. Theory. 30 (2013) 35-53. https://doi.org/10.1016/j.simpat.2012.07.013.

[22] C. VAN NGUYEN, S.K. SISTLA, S. VAN KEMPEN, N.A. GIANG, A. BEZOLD, C. BROECKMANN, F. LANGE, A comparative study of different sintering models for Al2O3, J. Ceram. Soc. Japan. 124 (2016) 301-312. https://www.jstage.jst.go.jp/article/jcersj2/124/4/124_15257/_article/-char/ja/.

[23] E.A. Olevsky, Theory of sintering: From discrete to continuum, Mater. Sci. Eng. R Reports. 23 (1998) 41-100. https://doi.org/10.1016/S0927-796X(98)00009-6.

[24] C. Manière, C. Harnois, S. Marinel, Porous stage assessment of pressure assisted
sintering modeling parameters: a ceramic identification method insensitive to final stage
grain growth disturbance, Acta Mater. 211 (2021).
https://doi.org/10.1016/j.actamat.2021.116899.

[25] B. Sarbandi, J. Besson, M. Boussuge, D. Ryckelynck, F. Barlat, Y.H. Moon, M.G. Lee, Anisotropic constitutive model and FE simulation of the sintering process of slip cast traditional porcelain, in: 2010: pp. 689–696. https://doi.org/10.1063/1.3457622.

[26] T. Kraft, H. Riedel, Numerical simulation of solid state sintering; model and application, J. Eur. Ceram. Soc. (2004). https://doi.org/10.1016/S0955-2219(03)00222-X.

[27] R.K. Bordia, R. Zuo, O. Guillon, S.M. Salamone, J. Rödel, Anisotropic constitutive laws for sintering bodies, Acta Mater. 54 (2006) 111–118. https://doi.org/10.1016/j.actamat.2005.08.025.

[28] O. Gillia, D. Bouvard, Phenomenological analysis of densification kinetics during sintering: Application to WC-Co mixture, Mater. Sci. Eng. A. (2000). https://doi.org/10.1016/S0921-5093(99)00621-8.

[29] T. Grippi, S. Béhar-Lafenetre, H. Friedrich, D. Haas, U. Schenderlein, S. Marinel, C. Manière, Grain growth modeling for gas pressure sintering of silicon nitride based ceramics, Mater. Today Commun. 34 (2023). https://doi.org/10.1016/j.mtcomm.2022.105189.

[30] T. Grippi, S. Béhar-Lafenêtre, H. Friedrich, D. Haas, U. Schenderlein, S. Marinel, C. Manière, Densification modeling for gas pressure sintered silicon nitride-based ceramics with Thermo-Optical dilatometry, J. Eur. Ceram. Soc. (2023). https://doi.org/10.1016/j.jeurceramsoc.2023.09.028.

[31] F. Raether, R. Springer, S. Beyer, Optical dilatometry for the control of microstructure development during sintering, Mater. Res. Innov. 4 (2001) 245–250. https://doi.org/10.1007/s100190000101.

[32] V.V. Skorohod, Theory of Nonlinearly Viscous and Plastic, 8 (1988) 621–626.

[33] C. Manière, E.A. Olevsky, Porosity dependence of powder compaction constitutive
parameters: Determination based on spark plasma sintering tests, Scr. Mater. 141 (2017) 62–66. https://doi.org/10.1016/j.scriptamat.2017.07.026.

[34] J. Besson, M. Abouaf, Rheology of Porous Alumina and Simulation of Hot Isostatic Pressing, J. Am. Ceram. Soc. (1992). https://doi.org/10.1111/j.1151-2916.1992.tb04479.x.

[35] C. Manière, L. Durand, A. Weibel, C. Estournès, Spark-plasma-sintering and finite element method: From the identification of the sintering parameters of a submicronic α-alumina powder to the development of complex shapes, Acta Mater. 102 (2016) 169–175. https://doi.org/10.1016/j.actamat.2015.09.003.

[36] Mohamed N. Rahaman, Sintering of Ceramics, in: Sinter. Ceram., CRC Press, 2007: pp. 1–43. https://doi.org/10.1201/b15869-2.

[37] M. Demartin, C. Hérard, C. Carry, J. Lemaître, Dedensification and anomalous grain growth during sintering of undoped barium titanate, J. Am. Ceram. Soc. 80 (1997) 1079–1084. https://doi.org/10.1111/j.1151-2916.1997.tb02949.x.

[38] C. Manière, G. Lee, E.A. Olevsky, Proportional integral derivative, modeling and ways of stabilization for the spark plasma sintering process, Results Phys. 7 (2017) 1494–1497. https://doi.org/10.1016/j.rinp.2017.04.020.

[39] Noboran, Graphite insulator properties, (n.d.). http://noboran.com/page.php?page=124/Graphite_Insulation__Boards.

[40] J.E.J.Ã. Staggs, H.N. Phylaktou, The effects of emissivity on the performance of steel in furnace tests, 43 (2008) 1–10. https://doi.org/10.1016/j.firesaf.2007.05.002.

[41] D. De Faoite, D.J. Browne, F.R. Chang-Díaz, K.T. Stanton, A review of the processing, composition, and temperature-dependent mechanical and thermal properties of dielectric technical ceramics, J. Mater. Sci. 47 (2012) 4211–4235. https://doi.org/10.1007/s10853-011-6140-1.

[42] J.J. Swab, A.A. Wereszczak, J. Tice, R. Caspe, R.H. Kraft, J.W. Adams, Mechanical and Thermal Properties of Advanced Ceramics for Gun Barrel Applications, Mater. Res.
(2005).

[43] K. Watari, K. Hirao, M.E. Brito, M. Toriyama, S. Kanzaki, Hot Isostatic Pressing to Increase Thermal Conductivity of Si 3 N 4 Ceramics, J. Mater. Res. 14 (1999) 1538–1541. https://doi.org/10.1557/JMR.1999.0206.

[44] A.F. Pan, W.J. Wang, X.S. Mei, K.D. Wang, W.Q. Zhao, T.Q. Li, Laser thermal effect on silicon nitride ceramic based on thermo-chemical reaction with temperature-dependent thermo-physical parameters, Appl. Surf. Sci. 375 (2016) 90–100. https://doi.org/10.1016/j.apsusc.2016.03.022.

[45] K. Watari, K. Sato, Y. Hotta, H. Nakano, K. Urabe, K. Ishizaki, S. Cao, K. Mori, Thermal conductivity of BeO-doped SiC ceramic at 5 to 1000 K, Key Eng. Mater. 264–268 (2004) 1083–1086. https://doi.org/10.4028/www.scientific.net/kem.264-268.1083.

[46] A. V. Blinder, A.S. Bolgar, V.Y. Petrovskii, Thermophysical properties of materials based on silicon nitride, Powder Metall. Met. Ceram. 34 (1996) 151–154. https://doi.org/10.1007/BF00559559.

[47] K. WATARI, Y. SEKI, K. ISHIZAKI, Temperature Dependence of Thermal Coefficients for HIPped Silicon Nitride, J. Ceram. Soc. Japan. 97 (1989) 174–181. https://doi.org/10.2109/jcersj.97.174.

[48] E. Olevsky, R. German, A. Upadhyaya, Effect of gravity on dimensional change during sintering—II. Shape distortion, Acta Mater. 48 (2000) 1167–1180. https://doi.org/10.1016/S1359-6454(99)00369-9.

[49] V.V. Skorohod, Rheological basis of the theory of sintering, Nauk. Dumka, Kiev. (1972).