Investigating the Ceramic Processes by Numerical Approaches

The considerable increase in computers performance has led in recent decades to a growing diffusion of numerical modelling as a valuable tool to support the manufacturing industry. Among the industrial sectors that may take advantage of this methodology there is the ceramic one. In fact, the production of large format ceramic slabs have introduced new technological challenges in a sector already subject to complex productdevelopment procedures, often based on a trial and error approach. For this reason, numerical simulations have become an actual and attractive method to overcome practical limitations and give a deeper insight into each phase of the production process. In this paper, advanced application of numerical modelling applied to the ceramic industry are presented, with special consideration for the production of ceramic slabs. Advantages and limitations of the adopted numerical techniques are discussed.


INTRODUCTION
The production process of ceramic slabs consists of many phases each characterized by its own variables and uncertainties. For instance, during the compaction phase the powder passes from a granular state to a fully dense one, thus making the mechanical properties of the granules and the characteristics of the tooling system the main factors on which the achievement of a defect-free and fully-dense green body depends. These aspects are, in turn, fundamental requirements to be reached prior to sintering, as dimensional changes and residual stresses in fired bodies are affected by the density distribution within the green parts as well as the way the cooling phase is performed [1].
In such context, numerical approaches can be used to support experimentations, by providing tools capable of predicting the behaviour of ceramic materials either during their production process [2][3][4] or during their designed operation [5]. The advantage of numerical approaches compared to traditional ones, which rely on trial and error procedures and on the technicians' expertise, consists in the possibility of virtually reproducing all the phases of the production process in order to verify their efficiency and reduce waste and time. In particular, numerical modelling allows to investigate each process involved in the production of ceramic materials and to identify its governing variables, thus improving the understanding of the processes themselves and paving the way to their optimization.
This study aims at showing the benefits in applying advanced numerical methods to the production process of ceramic slabs. The deposition, compaction, sintering and cooling phases are addressed and analysed. For each process, the main features incorporated into the mathematical model used for their simulation are discussed. Alongside the model description, the main variables involved in each phase will be identified and studied.

NUMERICAL STRATEGY AND RESULTS
In this section, the procedures for the numerical simulations of the deposition, compaction, sintering and cooling processes are presented. Numerical models have been developed to reproduce each of the abovementioned phases. In particular, the numerical approaches DEM (Discrete Element Method) and FEM (Finite Element Method), reliable to provide faithful results for all types of materials such as ceramics, composites and metals [5][6][7][8], have been used to reproduce the deposition process and the compaction and sintering processes, respectively. At the same time, experimental campaigns have been carried aiming at validating the obtained results.

Deposition Process
In order to simulate the deposition phase, the Discrete Element approach has been used. With this approach, powder particles are modelled as rigid spheres subjected to the gravity force and their deposition into a mould is dynamically simulated solving Newton's equations of motion [9]. The interactions between all the bodies (particles and boundaries) involved in the simulation is governed by contact laws, which depend on the material mechanical properties. Thus, the complex response at the macro-scale is obtained starting from the micromechanics driving the process at the scale of the grains. By using the DE Method provided by the opensource DEM platform YADE [10], it has been possible to investigate systematically different assemblies of spheres by varying their roughness and the Particle Size Distribution (PSD) (Figure 1

Compaction Process
Following the deposition, the spheres assembly is subjected to the compaction process. During this phase the granules bed undergoes large deformations accompanied by a significant densification, almost halving its original volume [11]. In order to develop a numerical model representative of this process, a different approach than DEM is required, which takes into account the spheres deformability. Thus, the geometry obtained from the DE model described above has been imported into a FE software [12]. The Finite Element approach, despite the higher computational cost if compared to DEM, gives greater flexibility as it takes into account complex physical phenomena between the grains occurring during the compaction phase. In particular, it assumes that the spheres are deformable elements whose behaviour can be determined given their mechanical properties. Thereby, the progressive deformation of the grains is taken into account and an estimation of the so-called "compaction curve" can be obtained. By using an elastic-perfectly plastic constitutive law, numerous numerical simulations have been performed on an assembly made up of around 1000 meshed particles placed into a cylindrical geometry ( Figure 2).
This modelling approach, based on FE simulations of the mechanical behaviour of an assembly of meshed particles in contact, is called Multi-Particle Finite Element Method (MPFEM). These simulations allowed to evaluate in a separate way the influence of some key parameters of the process, such as the friction between the grains and the yield stress of the material.
In parallel, to validate the proposed numerical model, experimental compaction tests in cylindrical moulds have been carried out ( Figure 3). By varying the PSD of the adopted powders and the dimensions of the mould, the influence of these variables on the compaction curve has been estimated. The compaction curves obtained from the numerical simulations show good agreement with the corresponding experimental curves. In addition, their comparison allowed the calibration of the MPFE-Model parameters in order to reproduce, with good accuracy, the process as a whole. To investigate further the evolution of the mechanical properties of the grains during compaction, a procedure based on [13] has been followed. An assembly of 50 spherical particles has been created using an algorithm implemented in YADE, which consists in generating randomly disposed particles in a cubic cell and growing their radii until contact between the particles is detected. At that point, contact forces develop at the spheres' interfaces as in the classical DEM procedure. Particles then rearrange until the stresses on the cell walls stabilise.
The obtained geometry has been imported into a FE code, each spherical particle has been meshed with approximately 3200 tetrahedral elements and the plates surrounding the spheres to form a cubic box have been modelled as rigid planes. The plates' side is 1 mm and the spheres radius is 0.13 mm. All spheres were modelled using a linear elastic constitutive law with Mises plasticity and an isotropic power-law strain hardening, with the material parameters taken equal to those reported in [13]. Contact parameters were also taken from the cited article. The so obtained assembly has been compacted under isotropic straining conditions (step-1) with subsequent unloading until low stress values (step-2).
Final relative densities ( ) of 0.90 and 0.98 have been reached at the end of this phase, to create a state of pre-consolidation within the assembly. Isotropic compaction has been obtained by shifting the rigid planes along their normal, thus assuming that stress and strain' principal directions remain coaxial through the simulations and allowing to study the response in the 3D principal stress space. After the first two steps, namely the strain-driven loading and unload-ding phases, a third step followed. In this phase, the sample has been reloaded under stress-controlled condition with a proportional loading path in a defined direction of the stress space to probe the material yield surface. Multiple loading directions have been tested by performing restarts analysis from the end of step-2. The onset of plasticity for each loading direction has been detected following the procedure explained in [12,13] which is based on the total dissipation within the numerical model. The obtained yield surfaces of the MPFE-Model for the two values of final relative density D are reported in Figure 4 in a deviatoric stress/mean stress diagram. The expression of the mean stress and deviatoric stress are defined by Equations 1 and 2.
( ) The explained procedure based on a meso-scale approach gives yield surface shapes that are comparable to experimental data available from the literature. It must be remarked that such models are characterised by a high computational cost as each restart analysis required around 10 hours to complete. A more efficient way to investigate the evolution of the material properties would be to use a continuum approach for the numerical representation of the powder material, strongly reducing the simulation time as well as the complexity of the model. Numerous studies such as [14] have been performed aimed at finding through the MPFEM the constitutive data of the material to be used in a continuum analysis. Among the continuum material models used for such analyses, the Gurson model [15] is one of the most popular. It assumes that the porous media is an aggregate of voids and a ductile matrix that is idealized as rigid-perfectly plastic and obeying the von-Mises yield criterion. The response of the media subjected to compaction depends on the matrix yield stress and porosity, and for a fully dense matrix the von-Mises yield criterion is recovered. Yield surfaces obtained with the Gurson model fitted to those given by the MPFEM simulations are shown in Figure 3. The shape of the yield surfaces is comparable, although an ad hoc hardening law should be included in the Gurson model to obtain a perfect agreement with the MPFEM results [14]. Finally, it must be noticed that a precise experimental validation of the relation existing between the grains and the macroscopic mechanical properties is particularly challenging in the ceramic sector, due to the difficulties encountered in the measurement of the grain properties. Nevertheless, such models can play a major role in improving our understanding of the compaction mechanism.

Sintering and Cooling
After compaction, the green body is subjected to sintering and lastly to cooling. Following these steps, the ceramic slab may present defects, often easily observable as cracks or deformations of the sides, or apparently invisible as the so-called residual stresses. Defects of the first type may be caused by a nonhomogeneous loading of the powders before compaction. In particular, it is observed experimentally that the density differences before sintering are usually much more pronounced than that after the sintering process completion. Such observation allowed to develop a simplified analysis conceived in order to reproduce the shape of the body at the end of the sintering process. Following such procedure, the body is treated as elastic and a space-varying volumetric shrinkage is imposed in order to lead all the body to the maximum densification at the end of the analysis. While such way of proceeding is clearly simplified and is not able to provide an evaluation of the stress distribution inside the slab (stresses are quickly dissipated thanks to the low viscosity of the material at high temperature), it is anyway able to provide a useful quantitative prediction of the deformed shape at the end of the sintering process.
In order to validate the proposed approach, an experimental campaign, aimed at verifying the effect of an inhomogeneous distribution of the powder during the deposition, has been performed. Experiments have been conducted using the experimental production line developed by SACMI, a leading company in the production of machines for the ceramic industry. Samples have been built by performing an inhomogeneous loading of rectangular tiles for which one of the quarters was compacted at a higher density with respect to the rest of the tile. Figure  4(a) shows one of such tiles (undeformed) while Figure  4(b) reports a comparison between the experimental and the numerical results, shown superimposed. The obtained numerical results show an excellent agreement with the experimental measurements in terms of deviation from the straightness of the sides. The deformations induced on the tile by the sintering and the remarkable overall tile shrinkage can be clearly observed as well. To what concerns the second type of defects, the attention has been focused on residual stresses. These are mainly caused by the different temperature occurring between the inner and outer part of the slab during cooling. In particular, at the beginning of the cooling phase the whole slab is at the temperature of maximum densification, which is around 1200 °C. As the temperature starts to decrease, the outer layers begin to cool down and thus contract, while the inside layers remain at a lower temperature.
Such phenomenon occurs when the material as a whole is still at very high temperature and extremely fluid, thus allowing the sliding between the outer layers and the innermost layers of the slab. At the end of the cooling process, the temperature is homogenized along the whole thickness inducing an "opposite" sliding compared to the previously described one. In this case, however, the phenomenon occurs at low temperatures (high viscosity), causing residual stresses to arise.
It must be noticed that the temperature difference between the outer layers of the slab and its inner part depends on both characteristics of the product (i.e. its thickness and material) and process parameters (i.e. the speed of the cooling process). If the stress induced by such temperature differences is too large and exceeds the material resistance when the material already assumed a brittle behaviour, cracks will occur. It is worth stressing that the cooling phase represents one of the most important bottleneck in the production process. It thus clearly appears that predicting stresses during cooling represents a very important aspect, which can be used in order to optimize the production process also taking into account the characteristics of the considered product. At the moment, such optimization is only seldom performed and it is based on the technicians' experience rather than on the base of numerical simulations.
To predict the magnitude of these stresses during cooling, a finite element model has been developed. First, the plate has been subjected to a uniform initial temperature of 1200°C, then the "virtual" cooling process has been started ( Figure 5). To quantify the effect of the speed at which the process takes place and following the procedure reported in [16], a parametric analysis has been carried out by applying two different cooling profiles, one fast and one slow, to the surfaces of the slab. In particular, (Figure 6) reports the evolution of the stress at the surface of the tile during the cooling process (fast cooling in this case). At the end of the cooling process, compression and tensile stresses arise on the external and internal surfaces respectively (Figure 7), with values increasing as the process speed increases, as expected.
It must be notice that, as the viscosity tends to infinity at the end of the cooling phase, when the material substantially assumes a purely elastic behaviour, the stress relaxation time tends to very high values so that they cannot be dissipated, so leading to residual stresses which might compromise the performance of the finite product. Such kind of defect is often observed as a reduction of the allowable stress in the final product with respect to the expected value and/or in difficulties in the cutting operations.
The model, currently in its validation phase, would allow to predict the magnitude of residual stresses given the material parameters and the cooling curves, so paving the way to a systematic optimization of the cooling phase.

CONCLUSION
The applications shown in this paper are just a few examples of how numerical simulation can be successfully applied to the ceramic industry, providing useful tools to support its rapid technological development. Despite the difficulties that often arise with the introduction of such techniques, especially in a complex sector as the ceramic one, the examples provided appear extremely encouraging and clearly show the great potentialities and advantages that these approaches can bring if combined with the currently used methodologies. In particular, it appears that, at the current stage, the optimization of the production process strongly relies on the technicians' expertise rather than on approaches based on simulation technologies. In the authors' opinion, numerous advantages might be obtained by the use of simulation technologies in the ceramic industry, following a path already well consolidated in many industrial sectors. A conspicuous amount of research still appears to be needed on such regard, in order to extend the use of numerical methods to the great complexity and multiphysics nature of the processes involved in the production of ceramic goods.