April 29, 2020

Numerical simulation of sintering for 3D-printed ceramics via SOVS model

Industry: Industrial Equipment | Software: pSeven | Company: Institute of Science and Technology

Introduction

Ceramic materials demonstrate excellent strength and hardness characteristics, high dimensional stability (low coefficient of thermal expansion), low density values, high resistance to abrasion and corrosion, and exceptional chemical stability. This unique combination of properties explains the growing demand for ceramic materials in different industrial fields including electronics, energy, mechanical and chemical engineering, and aerospace. However, the broad application of ceramic parts and components is restrained by complex and costly traditional manufacturing processes.

Today, rapid progress can be seen in ceramic 3D printing processes and their application to the manufacturing of complex-shaped parts. The major problem complicating the fabrication of high-quality parts is related to a substantial shrinkage of the parts at the sintering stage, resulting in significant distortions of the final part geometry compared to the expected result. It is a major impediment to the broad introduction of the technology. Hence, the application of mathematical modeling methods to simulate the behavior of a part during sintering has become a valuable tool to finely tune the process parameters.

Objective

The sintering models to describe changes in parts geometry using Finite Element Method are classified into phenomenological models and physical models. Skorohod-Olevsky Viscous Sintering (SOVS) phenomenological model has significant advantages as it can be easily implemented in numerical simulations and its most important input parameters can be determined by trial and error method, based on analysis of experimental data on sintering of test specimens.

The constitutive equation for the inelastic strain rate is expressed in SOVS as follows:

\( \dot{ \varepsilon }_{ij}^{in}=\frac{ \sigma ^{'}_{ij}}{2G_{p}}+\frac{ \left( \sigma _{m}- \sigma _{s} \right) }{3K_{p}} \delta _{ij} \) ,   (1)

where \( \sigma ^{'}_{ij} \) is the deviatoric stress, \( \sigma _{m}~ \) is the hydrostatic stress, \( \sigma _{s}~ \) is the sintering stress, \( \delta _{ij}~ \) is the Kronecker delta, and \( G_{p} \) and \( K_{p} \) are the viscous shear modulus and bulk modulus, respectively, which are defined in the following:

\( G_{p}= \left( 1- \theta \right) ^{2}AT^{n}\exp \left( \frac{B}{T} \right) \) ,  (2)

\( K_{p}=\frac{4 \left( 1- \theta \right) ^{3}}{3 \theta }AT^{n}exp⁡ \left( \frac{B}{T} \right) \) ,  (3)

where \( \theta = \left( 1- \rho \right) \) is the porosity factor and \( \rho \) is the relative density, \( A \) , \( B \) and \( n \) are material constants for the dynamic viscosity in a modified Arrhenius-type dependency for the fully dense ceramic material.

The sintering stress \( \sigma _{s} \) is expressed as follows:

\( \sigma _{s}= \frac{3 \gamma }{r} \left( 1- \theta \right) ^{2} \)     (4)

Where \( \gamma \) is the specific surface energy and \( r \) is the mean radius of ceramic powder particles \((d/2\)). For this study specific surface energy was considered as \( \gamma =0.9\frac{J}{m^{2}} \) , mean radius r was measured from 3D-printed Al₂O₃ specimens as 0.525 µm.

Thus, the objective of this use case was to identify unknown A, B and n parameters of SOVS model based on test data.

Solution

The specimens were 3D-scanned at every stage of the processing framework. The ceramic samples were thermally processed in accordance with the manufacturer’s suggested recipe, with variation implemented only in the part related to the sintering temperature and dwell values which were reduced to zero to evaluate the pure temperature effect. The following scenarios were tested with respect to sintering temperature: 1200–1700 °C with a step of 100 °C. After sintering, the specimens were scanned with the RangeVision 2M Pro 3D scanner to acquire a digital representation of the final geometry of the specimens and weighed with A&D HR250-AZG analytical scales. The porosity of every sample was estimated at both the beginning and end of the sintering process.

Figure 1. Evolution of Al₂O₃ material porosity after sintering.

The SOVS model was implemented in the ABAQUS modeling suite using the subroutine mechanism. 3D model of the part was created, just a quarter of the cylindrical test samples was considered for the calculations. Figure 2 displays the simulation results at a sintering temperature 1700 °C: Figure 2.a indicates the distribution of the magnitude of displacement in the sample at the end of the sintering process; Figure 2.b is the plot of the relative density evolution dynamics during the numerical simulations of the sintering process.

Figure 2. Results obtained at sintering temperature 1700 °C: (a) distribution of magnitude of displacements in sample at end of sintering process simulation and (b) plot of relative density evolution during numerical simulation of sintering process.

pSeven was used to identify the values of remaining parameters \(A\), \(B\) and \(n\) in SOVS model (2,3) based on experimental data for the specimens density. During the optimization cycle, pSeven performed a series of computations using ABAQUS/Standard solver to solve the direct numerical problem for A, B and n parameters. The pSeven workflow for the optimization problem is shown in Figure 3. A gradual modification and optimization of the parameters was made with the use of an internal pSeven algorithm that utilizes a combination of gradient-based methods for maximum acceleration of the convergence process. The absolute error was computed as the difference between the results of the current and previous iteration of the parameters; optimization continued until the error was less than 10^-4. The Root Mean Square Percentage Error (RMSPE) for difference between calculated and experimental density values was taken as the optimization criterion to be minimized, calculated as follows:

\( \delta =\min \sqrt[]{\frac{1}{6} \sum _{i=1}^{6} \left( \frac{ \rho _{i}^{\exp }- \rho _{i}^{sim}}{ \rho _{i}^{\exp }} \right) ^{2}} \cdot 100\% \),     (5)

where \(i\) – the number of the experimental point, \( \rho _{i}^{\exp } \) and \( \rho _{i}^{sim} \) – experimental and calculated values of relative density in the i-th point after sintering.

Figure 3. Solution scheme for optimization problem of finding values of A, B and n parameters for the dynamic viscosity.

The optimization strategy utilized the equation (1) for calculation of dynamic viscosity values, assuming that parameters A, B and n remain constant in the entire temperature range. At the first stage parameters A and B were varied while keeping parameter n constant and equal to 0.00, 1.00, 2.00, 3.00, 4.00, 5.00, consequently, for each start of optimization cycle. At the second stage all three pa-rameters were varied with the range of n between 4 and 5. The relation between optimization criterion and polynomial order n is shown in Figure 4(a). The best result with the minimum criterion of 0.64 % was achieved with 3 variable parameters (see Figure 5) within 187 iterations (Figure 4(b)). The best criterion values can be obtained with the following material parameters A = 1.5172 ・ 10^-10, B = 27716 and n = 4.26.

Figure 4. Convergence of optimization results: (a) Final optimization criterion as a function of n; (b) Convergence plot for three variable parameters A, B and n.

Figure 5. Comparison of calculated dependencies of final density on sintering temperature with ex-perimental values

Results

Table 1 shows the comparison between experimental shrinkages in the different directions and those predicted numerically. The analysis of results shows good convergence between calculated and experimental data. Especially good convergence can be observed at high temperatures at which parts are dwelled to achieve maximum relative density.

Sintering Temperature (°C)	Δx (%) (experimental)	Δy (%) (experimental)	Δz (%) (experimental)	Δ (%) (numerical)
1200	1.7±0.1	1.0±0.1	0.0*	0.4
1300	2.0±0.1	1.4±0.1	0.0*	1.1
1400	4.1±0.4	3.4±0.3	1.8±0.5	2.4
1500	6.4±0.2	5.8±0.2	5.8±0.7	4.9
1600	8.9±0.2	8.0±0.2	8.0±0.2	8.2
1700	11.9±0.4	11.0±0.3	12.8±1.7	11.6

It should be mentioned that the determined values of parameters A, B, and n are descriptive only for a particular tested class of ceramic feedstock. When modifications are introduced into the feedstock, i.e., powder fraction is changed or a novel feedstock is developed, a similar methodology must be applied to the new material to quickly estimate the suitable parameters for the description of the material densification.

The methodology demonstrated in this use case serves is an acceptable initial approach to studying the nature of the sintering process. In a simple form, it allows rapid scientific assessment of the physical effects of the sintering process in 3D-printed ceramics, even when large initial porosity is involved. The approach provides a walkthrough for estimating the parameters of the SOVS numerical model using pSeven optimization tools. The numerical model provides an acceptable fit to the experimental data, which enables the next stage of testing of its predictive power with more complex objects. From a technological perspective, the use of numerical simulations prior to the 3D printing of the ceramic object could significantly reduce the time and cost for the development and production of quality parts subject to strict tolerances. These tools could become a replacement for the methodology of using traditional shrinkage coefficients that is currently adopted by the majority of manufacturers.

Article is written based on the original publication "Numerical simulation of sintering for 3D-printed ceramics via SOVS model".