Skorohod-Olevsky Viscous Sintering Model Sensitivity to Temperature Distribution During the Sintering Process

This paper investigates the influence of temperature field non-uniformity on sintering simulation results using the Skorohod-Olevsky viscous sintering model.As a difference to previous studies, here a thermal transient analysis is performed to provide a detailed temperature field over the component within sintering time. Results obtained using uniform temperature distribution are compared to those obtained using a nonuniform distribution derived from a transient thermal analysis. Results are compared for different geometry sizes, that lead to different temperature non-uniformity levels. The study has shown that the temperature nonuniformity cannot always be neglected and should be considered as a possible source of modeling error.


INTRODUCTION
Sintering represents a phase in component production in which a sample produced by powder compaction (socalled green body) is subjected to combined high temperature and pressure, in order to obtain the final product. During the process, cohesive forces are developed and porosity is gradually reduced, or completely eliminated, with consequent change in sample's dimensions. To keep the shape and dimensions of the final sample under control, its local shrinkage, and potential geometry distortions during the sintering process should be predictable. The knowledge on the subsequent shrinkage should be employed when designing a corresponding green body. In order to avoid a trial and error approach in production, it is desirable to have a reliable numerical model to simulate the whole process.
To these purposes numerical simulations based on finite element method are frequently employed both for the initial powder compaction phase and for the subsequent sintering phase [1,2]. In general, two different approaches are used to model mechanical response of the material: (1) micro-mechanical modeling which is apt to capture the phenomena that are taking place at the lower scale, and (2) continuum based modeling, by relying on phenomenological constitutive models [3,4]. In the past two decades, the latter approach has gained popularity, since being more appropriate form practical engineering purposes, as the problem can be solved over one scale only.
The main feature of all continuum based (i.e. phenomenological) models is that the shrinkage duringpowder compaction and subsequent sintering is modeled by inelastic strain. A popular phenomenological model that is used for sintering simulations is the Skorohod-Olevsky viscous sintering (SOVS) model. It is a visco-plastic model based on rheological sintering theory [5,6], modified for the use within continuum mechanics context [7]. The original model and its subsequent extensions and improvements have been successfully used for modeling of the sintering process [8,9].
The model is primarily intended to be used within a loosely coupled thermo-mechanical analysis using temperature as an auxiliary variable with a precalculated evolution (e.g., uniform temperature change or a result of a thermal analysis). A lot of work has been done on the application of the model for different materials, implementing diverse numerical solution algorithms and evaluating the influence of spatial and temporal discretization on the analysis convergence performance and the error estimations [10]. However, most of the published papers involving sintering simulations using the SOVS model do not deal with temperature field in detail, also in the case of multilayers sintering [11,12], where it is expected to have stronger temperature nonuniformity. Commonly, a uniform temperature evolution of the sintered body is defined based on sintering time, temperature and heating and cooling rates.This temperature profile isadopted as an input for sintering simulation.
In this work, the sensitivity of SOVS model with respect to temperature distribution in the sintered part is studied. Uniform temperature distribution is compared to a more realistic case considering distribution resulting from a detailed thermal analysis. Depending on the heating method, component size and material properties, temperature distribution in the sample during sintering can range from nearly uniform to significantly nonuniform. The work presented here attempts to capture the influence of these differences to the final result of sintering simulation.In what follows first a brief description of the adopted constitutive model is presented, in tensor notation.

CONSITUTIVE MODEL
Typically used approach in plasticity models, employed also for the SOVS model consist in additive decomposition of the total strain into its elastic and inelastic part. Written in a rate form it reads: The rates here refer to fictitious time, hence should be viewed as the small changes in strain value, rather than time derivatives. This formulation is typical for plasticity (see e.g. [13]) and is adopted also here for the inelastic behavior.
For the elastic part a linear elasticity is assumed, which is consistent with the application of sintering simulations where only small elastic deformation are encountered. The stress is hence related to the elastic part of the strain using Hooke's law: With μ and λ being Lame's parameters (i.e. elastic constants of the material) and δ ij is the Kronecker delta.For consistency in writing equation (2) is expressed in rate form. Note that for linear elastic behavior the rate form of equation matches the one expressed in the absolute value of strain. Hence the elastic constants in (2) are related to more commonly used Young's modulus and Poisson's ratio through: ( ) ( )( ) During the sintering phase the relative density is continuously changing giving rise to the inelastic deformation. Therefore within this phenomenological model the change in relative density is modeled through the inelastic strain, or specifically, through its volumetric part. Thus, the relative density is defined as the ratio of current material density and theoretical maximum density: The complete inelastic strain rate within the model is governed by the following constitutive equation: where σ s is the effective sintering stress (Laplace pressure), G is the effective shear viscosity and K is the effective bulk viscosity. Two auxiliary state variables in material point are also used:temperature T and relative density ρ. While the latter is computed on the basis of the inelastic strain, through equation (3), the former is either pre-defined, or, like in this study, results from the thermal analysis of the component. Constitutive parameters entering into equation (5) (i.e. the effective viscosities G and K ) depend on relative density ρ through normalized shear and bulkviscosities φ and ψ, but also depend on temperature through shear viscosity of the fully dense material : Most commonly, these functions are defined as: where a i , b i and c i are model parameters, whilelocal sintering stress σ s0 is given by: where α is thesurface tension and r the average grainsize.
There have been several efforts to improve the model by developing new forms of these functions. Majority of them have been aimed at improving model response for certain materials by taking into account some specific material features such as inelastic deformation mechanism or grain growth [14]. One notable improvement has been made by using an Arrhenius-type shear viscosity function [8], given by the following relation: This is the form of the shear viscosity function used in this work, and the corresponding model will be referred to as the modified SOVS (mSOVS) model.
It is worth noting that although this model uses a variable density, the model itself is intended to be used within quasi-static simulations.Hence, no inertial forces are computed within the calculations. Relative density is used as an auxiliary state variable, which is computed purely on the basis of inelastic part of the strain, as previously outlined. Therefore the model represents a pure phenomenological visco-plastic mechanical model.

METHODOLOGY
A common mechanical component shape is considered (Figure 1) as a case study for performing sintering simulations within this work. Only one dimension is given in the figure to provide informative basis of the global size of the component.Sintering was simulated by employing the open source finite element software Code Aster [15]. As the above described constitutive model does not exist as already available feature within the employed software, for the purposes of this study it was implemented as User subroutine. A modified version (mSOVS) was used, while the numerical imple-mentation was done in accordance to details specified in work [16]. It is the most widely used numerical implementation of the considered model. Within this study a fully implicit version was implemented. For details the reader should refer to [16].

Figure 1. Sample geometry
Sintering temperature profile considered within the simulation started from ambient temperature of 20°C, applying first a heating with the rate of 50°C/min until reaching sintering temperature of 1120°C. This temperature was further kept constant for the duration of 2 hours, followed by a subsequent cooling with the same rate of -50°C/min until reaching again the ambient temperature of 20°C. Distribution of relative density at the beginning of sintering was assumed to be changing linearly over the height with average value of 0.8. Alternatively, distribution resulting from powder pressing simulation can be used (e.g. like the one provided from simulations with modified Drucker-Prager cap model -see [2]). Such an approach clearly represents a more realistic scenario. However, in this work the focus was on the influence of non-uniform temperature distribution on the final deformation field over the sintered component, so the other causes were excluded from the consideration. To determine the model sensitivity to temperature distribution, two groups of simulations have been performed: • sintering simulation using a uniform temperature field over the whole specimen as an input (Figure 2), • sintering simulation using a non-uniform temperature field resulting from a transient thermal analysis as an input. Result of transient thermal analysis clearly depends on the shape of the component considered but also on its size, as for the larger component if would be generally required a longer time to reach steady temperature distribution. It is therefore expected to have a stronger nonuniformity of the temperature field for larger scale components. Such influence is also investigated within this study as the additional comparative basis by performing simulations considering the uniformly scaled sample geometry. Specifically five different scales were considered: 2, 4, 8 and 16 times scaled the reference one visualized in Figure 1.
In the transient thermal analysis, a heat transfer boundary condition has been applied to all outer surfaces of the sample through the following equation: where λ is material thermal conductivity, n is surface normal, h is heat transfer coefficient and T ext is sintering profile temperature according to Figure 2. Heat transfer coefficient has been chosen so that the resulting temperature field for the reference size sample (i.e. with dimensions visualized in Figure 1) is nearly uniform, but is significantly non-uniform for the largest sample considered. Its value is assumed to be equal to: The following mechanical and thermal properties of the material have been used, according to the values that can be found in the literature for alumina (Al 2 O 3 ): • linear thermal expansion coefficient As the adopted model is fully phenomenological one, the material is considered as continuum, so the alumina is not treated as a powder.
An initial relative density distribution, characteristic for die pressed green bodies, shown in Figure 3has been taken for all simulations the same.

Figure 3. Initial relative density distribution for all simulations
In order to reduce computing time both symmetries were exploited so only one quarter of the sample was modeled (see Figure 4). Symmetry boundary conditions were imposed on cutting planes, without additional constrains to simulate a free standing sample. The same finite element mesh was used both for thermal and mechanical simulations with 15948 tetrahedron finite elements.

Transient thermal analysis results
Thermal analyses revealed that for the nominal geometry (i.e. the one visualized in Figure 1), the temperature distribution is nearly uniform over the sample, throughout the whole sintering time. Similar result is obtained also when the sample is scaled 2 times. The difference starts to be noticeable for the size scaled by 4 times. For this sample the temperature distribution resulted in internal areas with about lower temperatures than nominal according to the sintering profile. Figure 5 visualizes the distribution of temperature over the sample after 19minutes from the beginning of heating.
To assess the overall difference in temperature uniformity for diverse sample sizes, an average sample temperature with respect to simulation time for different model scales is given on Figure 6.The average temperature is calculated as a mean value of nodal temperature values. Averaging by considering weight coef-ficients according to element volumes was not considered here due to relatively uniform element sizes within the adopted mesh(see Figure 4).The figure evidences that for the largest scale considered the uniform temperature distribution is reached only at the very end of sintering time. Such temperature profile will lead to earlier development of cohesive forces for the areas closer to the outer surfaces, followed by nonunique shrinking over the specimen, potentially giving rise to residual stresses and shape distortion. It is expected that these differences will eventually lead to more pounced residual deformations upon finishing the sintering process. Clearly in such scenario the simulation with ideally uniform temperature distribution may yield unrealistic results.

Sintering simulation results
The change in the average sample's relative density during the sintering process for different model scales is given on Figure 7. As expected, larger samples took longer time to reach sintering temperature which resulted in lower final density, i.e. higher final porosity. Average final density error induced by approximating the temperature distribution as uniform with respect to model scale is shown on Figure 8.
In the following, simulation result difference induced by approximating temperature field as uniform will be considered as an error. For relatively smaller sample sizes, the final sample average density error is negligible . It is important to note that for larger scale samples larger values of relative density (i.e. smaller amount of porosity) is reached closer to the surface, so the averaged relative density over the whole sample is not providing detailed enough insight into the final quality of the component. Figure 9 visualizes the distribution of relative density over one larger component of here considered (i.e. the one scaled by 8 times the nominal one). The result visualized in the figure shows that local relative density over the sample varies from 0.84 to 0.94, evidencing the areas with larger level of porosity (i.e. lower value of relative density).Such result provides an important indication that the average relative density is not sufficient to judge upon achieved quality of sintering, as there can be significant difference between local and average relative density. Here specifically, the smallest value is about 0.84, while the average relative density is equal to about 0.91 (see Figure 7). The difference is even bigger for the largest scale component, not shown here for brevity. Over the component it is noticeable that a general gradient of reduction of relative density with height is still preserved to a certain extent, which is due to the difference in starting relative density (see Figure 3).
Another important aspect of quality verificationfor the final component that can be assessed through the simulation of sintering concerns controlling of the final shape and dimensions ( Figure 10). Non-uniform temperature distribution during sintering in general leads toearlier densification of some areas and different shape distortion when compared to the case of uniform temperature distribution ( Figure 11).
Average final displacement error induced by approximating the temperature distribution as uniform with respect to model scale is shown on Figure 12.

CONCLUSION
In this work numerical simulations of sintering process were conducted. Combined thermal and mechanical simulations of the sintering process were performed by employing a phenomenological constitutive model. This numerical setup provided a framework to study the influence of temperature non-uniformity during the sintering, as opposed to commonly used simplification of uniform temperature distribution over the sample. The study revealed a general conclusion that this non-uniformity becomes rather important for larger scale components, particularly if one is interested in local distribution of porosity and not just the averaged relative density. The obtained results also corroborate to some more specific conclusions which are outlined in the following.
It was evidenced that the final sample's average density is not affected much and for smaller sized components can be even neglected. However by looking to this criterion only, a misleading conclusion can be obtained. The results of here conducted numerical simulations revealed that a slightly different relative density distribution due to a nonuniform temperature field during sintering can influence final sample shape (Figures 10 and 11). This can become rather important in design of mechanical parts with tighter tolerances, as the residual deformation can lead to violation of these acceptable dimension ranges. Shape distortion can be captured well by the numerical simulations employed here, while the result will be more realistic if the thermal analysis, apt to provide the temperature distribution over the whole sintering time, is also considered.
Another important aspect representslocal residual stresses and local porosity distribution. It can be stated that in general, larger scale components provoke less uniform distribution of porosity over the sample. The same can be claimed for the shape distortion, which in turn provokes a larger self-equilibrated stresses in the structure. Stresses are varying locally in the component and their magnitude cannot be captured well just by observing a global quantity, averaged over the sample, like relative density. The knowledge of these stresses may be important as they could lead to premature failure of the component due to fatigue. If a local distribution of any state quantity is in focus, the thermal analysis prior to sintering simulation becomes rather important, as it provides more realistic temperature distribution over the sample.
Temperature field non-uniformity depends on many parameters, and not only on sample's size and geometry. As an example heating method, heating rate and material properties will all influence the resulting temperature field. Transient thermal analysis used here is a computationally inexpensive way of determining average sample temperature during sintering, by considering all the factors. Deviation of calculated average sample temperature in comparison with the chosen sintering profile is a good starting point for assessing possible error induced by approximating temperature distribution as a uniform one.

ACKNOWLEDGMENT
Results presented in this paper are related to the project that received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 955944 -RE-FRACTURE2.