Charles Manière
email: charles.maniere@ensicaen.fr

Christelle Harnois

Sylvain Marinel

3D printing of porcelain: finite element simulation of anisotropic sintering

Keywords: 3D printing, Finite Element Simulation, Anisotropy, Mesoscale Simulation, Sintering

Sintering is a key step that determines the mechanical performances and shape quality of 3D printed ceramics. In this domain, there is an increasing interest in complex thin structures and low fill density shapes. However, these complex structures may suffer from distortions by gravity forces developed on the light structure, friction with supports, shrinkage anisotropy or thermal gradients. Finite element simulation of the sintering process is then of high importance for predicting the specimens distortions. This paper focuses on the identification of the model parameters in the special case of anisotropic sintering of porcelain with final stage swelling. This study points out the printed specimen high resistance to high temperature shear deformation due to the alignment of the printed rods. A mesoscale simulation has been carried out to explain this mechanism. The resulting model was tested for the sintering of a thin wall cup shape.

Introduction

Ceramic 3D printing opens an interesting window for the design of high-resolution materials [START_REF] Chartier | Fabrication of Millimeter Wave Components Via Ceramic Stereo-and Microstereolithography Processes[END_REF][START_REF] Deubel | Direct laser writing of threedimensional photonic-crystal templates for telecommunications[END_REF] with complex structures based on multiscale/bioinspired shapes [START_REF] Li | Bioinspired Functional Surfaces Enabled by Multiscale Stereolithography[END_REF][START_REF] Yang | 3D-Printed Biomimetic Super-Hydrophobic Structure for Microdroplet Manipulation and Oil/Water Separation[END_REF], low density/lightweight structures [START_REF] He | Support-Free Ceramic Stereolithography of Complex Overhanging Structures Based on an Elasto-viscoplastic Suspension Feedstock[END_REF][START_REF] Li | Additive manufacturing of lightweight and highstrength polymer-derived SiOC ceramics[END_REF] and advanced composite structures with improved mechanical/functional properties [START_REF] Zhang | Stereolithography-based additive manufacturing of lightweight and high-strength Cf/SiC ceramics[END_REF][START_REF] Feilden | 3D Printing Bioinspired Ceramic Composites[END_REF]. Ceramic additive manufacturing involves different methods with resolution from 0.4-1 mm for robocasting [START_REF] Peng | Ceramic Robocasting: Recent Achievements, Potential, and Future Developments[END_REF] (or direct write printing) to 10-100 µm for stereolithography [START_REF] Chartier | Rapid Prototyping of Ceramics[END_REF]. These processes produce green specimens by extrusion and ceramic resin photopolymerization, respectively [START_REF] Chen | 3D printing of ceramics: A review[END_REF][START_REF] Deckers | Additive Manufacturing of Ceramics: A Review[END_REF][START_REF] Wang | Review of additive manufacturing methods for high-performance ceramic materials[END_REF]. 3D printed green specimens are often made of low-density highly complex shapes and a complex support material network [START_REF] Jiang | Support Structures for Additive Manufacturing: A Review[END_REF]. The main motivation is to obtain lightweight materials with an optimal mechanical resistance suitable for a target application [START_REF] Du | 3D printing of ceramic-based scaffolds for bone tissue engineering: an overview[END_REF]. However, these complex shapes are fragile all along the different processing steps, i.e. the printing, debinding and sintering.

During the printing, thin shapes are very sensitive to gravity slumping. This phenomenon is controlled by the ceramic paste rheology and the printing conditions [START_REF] Miyanaji | Process development for green part printing using binder jetting additive manufacturing[END_REF][START_REF] Barki | Linking Rheology and Printability for Dense and Strong Ceramics by Direct Ink Writing[END_REF]. The debinding is also a key step where cracks and swelling may appear if the temperature cycle is not carefully adapted to the organic phase decomposition [START_REF] Bae | Influence of Residual Monomer on Cracking in Ceramics Fabricated by Stereolithography[END_REF][START_REF] Li | Dental ceramic prostheses by stereolithography-based additive manufacturing: potentials and challenges[END_REF]. Finally, thin shapes may suffer from anisotropy [START_REF] Manière | Modeling sintering anisotropy in ceramic stereolithography of silica[END_REF][START_REF] Tarabeux | Development of a numerical simulation model for predicting the curing of ceramic systems in the stereolithography process[END_REF] and significant distortions [START_REF] Withell | Porous ceramic filters through 3D printing[END_REF][START_REF] Ohji | Additive manufacturing of ceramic components[END_REF] at the sintering stage due to various causes such as gravity, friction with the support [START_REF] Alvarado-Contreras | Modeling of gravity-induced shape distortions during sintering of cylindrical specimens[END_REF], high material viscosity for liquid phase sintering [START_REF] Alvarado-Contreras | A continuum approach for modeling gravitational effects on grain settling and shape distortion during liquid phase sintering of tungsten heavy alloys[END_REF][START_REF] Shanshan | An Experimental Study of Ceramic Dental Porcelain Materials Using A 3D Print (3DP) Process[END_REF], swelling (or bloating) at the final stage [START_REF] Riley | Relation of Chemical Properties to the Bloating of Clays[END_REF][START_REF] Manière | Swelling negation during sintering of sterling silver: An experimental and theoretical approach[END_REF] and thermal gradients [START_REF] Manière | Fluid dynamics thermo-mechanical simulation of sintering: Uniformity of temperature and density distributions[END_REF].

Cantilever shapes without support materials are also highly sensitive to shape distortions during the sintering [START_REF] Lynnora | Mitigating Distortion during Sintering of Binder Jet Printed Ceramics[END_REF].

For selective laser melting/sintering (SLM/SLS technology), the multiphysics finite element tool is very precious to evaluate the mesoscale and macroscopic behavior of complex structures mechanical properties [START_REF] Eshraghi | Micromechanical finite-element modeling and experimental characterization of the compressive mechanical properties of polycaprolactonehydroxyapatite composite scaffolds prepared by selective laser sintering for bone tissue engineering[END_REF] including residual stresses, melting/sintering phenomena [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramic materials[END_REF] and shrinkage issues [START_REF] Gajera | Experimental investigation and optimization of direct metal laser sintering process for shrinkage rate using CL50WS material[END_REF][START_REF] Fotovvati | Size effects on geometrical accuracy for additive manufacturing of Ti-6Al-4V ELI parts[END_REF]. The finite element simulation of the sintering process is highly important to predict the distortion resistance of the different complex shapes studied by 3D printing. In particular, such simulation can predict the mechanical resistance of the green specimen during its sintering, allowing adjusting the different strategies of filling density and/or the structure of the sacrificial support materials. This work is dedicated to the sintering of 3D printed porcelain obtained by robocasting [START_REF] Yu | Extrusion-based additive manufacturing of yttriapartially-stabilized zirconia ceramics[END_REF]. Liquid phase sintering of porcelain is a very sensitive process which requires a careful control of the cycle, in particular at the end of the sintering stage, where swelling may occur due to trapped gas in pores generating a force opposite to the capillarity sintering driving force [START_REF] Sánchez | Porcelain tile: Almost 30 years of steady scientific-technological evolution[END_REF]. This phenomenon is particularly important if the final stage temperature is too high [START_REF] Martín-Márquez | Effect of firing temperature on sintering of porcelain stoneware tiles[END_REF]. However, the main phenomenon studied in the present paper is the sintering shrinkage anisotropy resulting from the extruded rods stacking strategy, printing conditions and consequently from the porosity generated [START_REF] Peng | Ceramic Robocasting: Recent Achievements, Potential, and Future Developments[END_REF][START_REF] Chartier | Manufacturing of ceramic parts by additive manufacturing technologies[END_REF]. The impact of printing architecture, shape scale and sintering conditions on anisotropy will be investigated. The previously developed method for the identification of anisotropic behavior will be used to determine the sintering model anisotropic parameters [START_REF] Manière | Modeling sintering anisotropy in ceramic stereolithography of silica[END_REF]. Thereafter, the resulting model will be used in a finite element code to calibrate the viscous behavior of the printed porcelain specimen via a deflection test and a mesoscale simulation. Finally, the simulation will be used to predict the anisotropic sintering shrinkage behavior of printed complex shapes.

Theory and calculations

The sintering model is based on the continuum theory of sintering [START_REF] Olevsky | Theory of sintering: from discrete to continuum[END_REF]. In the following, the sintering general equations will be defined first; then, the modeling case of sintering anisotropy in the presence of porosity will be described.

Sintering continuum model

The continuum theory of sintering defines the sintering behavior of a continuum compressible medium by the following stress and strain rate tensor equation [START_REF] Olevsky | Theory of sintering: from discrete to continuum[END_REF][START_REF] Bordia | Current understanding and future research directions at the onset of the next century of sintering science and technology[END_REF].

𝜎 = 𝜎 𝑒𝑞 𝜀̇𝑒 𝑞 (𝜑𝜀̇+ (𝜓 - 1 3 𝜑) 𝑒i) + (𝑃 𝑙 -𝑃 𝑠 )𝕚 (1)
The right term 𝑃 𝑙 𝕚 is the capillarity force resulting from sintering stress that can be expressed theoretically via Skorohod equation (2) [START_REF] Skorohod | Rheological basis of the theory of sintering[END_REF].

𝑃 𝑙 = 3𝛼 𝑟 (1 -𝜃) 2 (2).
In equation ( 1), the equivalent stress 𝜎 𝑒𝑞 and strain rate 𝜀ė 𝑞 ratio corresponds to the equivalent dense phase behavior. As pressureless sintering has low sintering stresses, the dominant sintering matter transport mechanisms are assumed to be diffusive or viscous mechanisms.

These mechanisms are linear and obey linear viscous equation (3) [42].

𝜎 𝑒𝑞 = 2𝜂𝜀̇𝑒 𝑞 (3).
For porcelain sintering which is based on liquid phase sintering (also called vitrification sintering), the mechanism is typically linear and also obeys equation [START_REF] Li | Bioinspired Functional Surfaces Enabled by Multiscale Stereolithography[END_REF].

The temperature dependence of the material viscosity 𝜂 is typically [START_REF] German | Sintering Theory and Practice[END_REF]:

2𝜂 = 𝜂 0 𝑇exp ( 𝑄 𝑅𝑇 ) (4).
Equation (1) general form for pressureless sintering is then:

𝜎 = 2𝜂 (𝜑𝜀̇+ (𝜓 - 1 3 𝜑) 𝑒i) + (𝑃 𝑙 -𝑃 𝑠 )𝕚 (5)
.

The porosity function 𝜑 and 𝜓 (shear and bulk modulus respectively) can also be theoretically defined via Skorohod's equations [START_REF] Olevsky | Theory of sintering: from discrete to continuum[END_REF][START_REF] Skorohod | Rheological basis of the theory of sintering[END_REF]:

𝜑 = (1 -𝜃) 2 (6)
𝜓 = 2 3
(1-𝜃) 3 𝜃 [START_REF] Zhang | Stereolithography-based additive manufacturing of lightweight and high-strength Cf/SiC ceramics[END_REF].

In order to link the volume change rate (𝜀ẋ + 𝜀ẏ + 𝜀ż) and the porosity elimination rate 𝜃 , the mass conservation equation is employed.

Implementing equations (5,6,7 and 9) in a finite element code allows simulating the pressureless sintering. The main sintering parameters to determine are the following 𝜂 0 /𝛼, 𝑄 and the initial particle radius (see details in [START_REF] Manière | Modeling sintering anisotropy in ceramic stereolithography of silica[END_REF]). The shear and bulk moduli may also be determined experimentally [START_REF] Manière | Porosity dependence of powder compaction constitutive parameters: Determination based on spark plasma sintering tests[END_REF], however using Skorohod equations [START_REF] Skorohod | Rheological basis of the theory of sintering[END_REF], porous behavior can be approximated reducing drastically the identification methodology.

𝜎 = 2𝜂 (𝜑𝜀̇+ (𝜓 - 1 3 𝜑) 𝑒i) + (𝑃 𝑙 -𝑃 𝑠 )𝕚 (10)
.

Anisotropic pressureless sintering model

The 3D printing process implies a preferential interlayer porosity distribution [START_REF] Manière | Modeling sintering anisotropy in ceramic stereolithography of silica[END_REF][START_REF] Tarabeux | Development of a numerical simulation model for predicting the curing of ceramic systems in the stereolithography process[END_REF].

Consequently, the sintering shrinkage of the parts is rarely isotropic and may experience an anisotropic behavior. This anisotropic behavior can be inserted at the level of the viscosity 𝜂, shear and bulk moduli ( 𝜓and𝜑) or sintering stress 𝑃 𝑙 . However, in opposition to the nonlinear sintering law, the linear sintering equation does not distinguish the two latter terms that are generally defined together as shear and bulk viscosity [START_REF] Manière | Modeling sintering anisotropy in ceramic stereolithography of silica[END_REF][START_REF] Kraft | Numerical simulation of solid state sintering; model and application[END_REF]. Therefore, the anisotropy behavior can be inserted at the level of the viscosity (𝜂 𝑧 , 𝜂 𝑟 ) having then a different value in the layers plane (R axis) and in the normal layers direction (Z axis).

2𝜂 𝑟 = 𝜂 0𝑟 𝑇exp ( 𝑄 𝑅𝑇 ) (11)
.

2𝜂 𝑧 = 𝜂 0𝑧 𝑇exp ( 𝑄 𝑅𝑇 ) (12)
.

The anisotropic model parameters can be identified by adjusting the simulated sintering response to the dilatometry specimen diameter and height to obtain: 𝜂 0𝑟 , 𝜂 0𝑧 and 𝑄. The finite element simulation code employed was Comsol Multiphysics.

Experiment and method

This study is divided in three main stages including (i) the modeling parameters identification from dilatometry (ii) the correction of the model based on deflection sintering tests combined with the simulation tool and (iii) the final comparison of the resulting modeling approach on a complex shape.

Samples preparation

The sintering dilatometry tests were conducted on the apparatus "Setaram TMA92" with 3D printed porcelain cut in 6.7 mm height and 5.7 mm diameter specimens. The robocasting 3D printing method was based on "delta wasp 2040 clay". An aqueous porcelain ceramic paste consisting of "CERADEL PT010B" with 5 % of water addition was used for the robocasting process with the following printing conditions: print speed 50mm/s, nozzle size 1.2 mm, extrusion gas pressure 4 bar, layer height 0.5 mm, bottom/top thickness 2 mm, shell thickness 2.4 mm, filling density 100 %. In a previous article [START_REF] Manière | 3D printing of porcelain and finite element simulation of sintering affected by final stage pore gas pressure[END_REF], a dilatometric study was carried out showing the swelling phenomenon and allowing to determine an optimum sintering cycle providing high density while minimizing the pore gas pressure swelling phenomenon. This cycle is 1K/min to 1200°C with 2h holding. This thermal cycle has been used for all the experiments of this study, including for the dilatometric measurements. In order to determine the anisotropic sintering shrinkage magnitude for different experimental parameters, three runs of tests were done to investigate the effect of sintering heating rates, volume and filling ratio. The identification of the sintering parameters using equations ( 11and 12) will be done by the adjustment of the specimen dimensions.

Bar deflection tests and anisotropic simulation

In order to test and adjust the sintering parameters identified by dilatometry, the numerical tool is first compared to a bar deflection test. The latter dimensions are 6 × 6 × 60 mm. The bar is placed on the top of an alumina cup with a 45 mm inner diameter for allowing the specimen slumping during sintering. The configuration is presented in figure 1a. Figure 1b illustrates the distortion after the sintering (1K/min to 1200°C, 2h holding). A deflection of 0.5 mm is observed. This experimental deflection will be compared to the finite element simulation result.

Sintering anisotropy characterization

Figure 2 shows the results of the post-sintering anisotropy study. The anisotropy is characterized by measuring the ratio of the final and initial height and diameter. If the sintering is isotropic, the ratio hf/hi and df/di must be the same, and different if anisotropic shrinkage occurs. In a previous study [START_REF] Manière | 3D printing of porcelain and finite element simulation of sintering affected by final stage pore gas pressure[END_REF], sintering with different heating rates has been explored by dilatometry on the "10 mm" specimens. The results presented in figure 2a, show an anisotropic shrinkage behavior with hf/hi close to 0.83 and df/di close to 0.9. The heating rate does not have any noticeable impact on these ratios. Very similar results are obtained for the scalability and filling ratio tests illustrated in figure 2b and figure 2c respectively. Bigger specimens and high void ratios in the volume of the specimen do not generate significant differences in the anisotropy behavior. In overall, the average hf/hi ratio is 0.89 ± 0.01 and df/di ratio is 0.82 ± 0.01. This reproducibility of the anisotropic sintering shrinkage indicates the latter must be originated from the layered printing conditions rather than the heating cycle, the specimen size or its architecture (filling ratio). This means that the sintering shrinkage can be predicted by the finite element simulation after a careful identification of the anisotropic sintering conditions. In below section the identification process begins with the dilatometry test. reproduce these anisotropic sintering responses. As a starting point, the isotropic model parameters employed in ref [START_REF] Manière | 3D printing of porcelain and finite element simulation of sintering affected by final stage pore gas pressure[END_REF] will be used and modified to reproduce the anisotropic behavior. Two modeling approaches are investigated; one will consider the trapped gas pressure in the closed pores at the final stage [START_REF] Manière | 3D printing of porcelain and finite element simulation of sintering affected by final stage pore gas pressure[END_REF] and will be referred as "pore gas pressure model". The other modeling approach will consist in the use of a threshold law on the bulk modulus [START_REF] Chen | 3D printing of ceramics: A review[END_REF] to explain the sintering slowdown at the final stage (by a critical porosity 𝜃 𝑐 ) instead of considering the trapped gas pressure in pores. This approach will be referred as "threshold model".

𝜓 = 2 3
(1-𝜃) 3 (𝜃-𝜃 𝑐 ) [START_REF] Chen | 3D printing of ceramics: A review[END_REF].

The two sintering models are reported in figure 3 after the adjustment of their parameters.

Both succeed to reproduce the sintering response of the printed sample and have in common the activation energy of 300.3 kJ/mol and 𝜂 0𝑧 = 3.00𝐸 -3𝑃𝑎. 𝑠 . Then, the "pore gas pressure model" considers the function 𝑃 𝑠 determined in ref [START_REF] Manière | 3D printing of porcelain and finite element simulation of sintering affected by final stage pore gas pressure[END_REF], the shear modulus expression [START_REF] Zhang | Stereolithography-based additive manufacturing of lightweight and high-strength Cf/SiC ceramics[END_REF] and:

𝜂 0𝑟 = 𝜂 0𝑧 + 0.001𝜃 2.5 (12)
.

The porosity function in the expression of 𝜂 0𝑟 is a variable tending to zero at full density. The latter function gradually suppresses the anisotropy behavior as the porosity decreases. This function considers the porosity distribution (higher interface porosity) as the origin of the anisotropic behavior. Consequently, at full density, the anisotropic model tends to a pure isotropic creep model [START_REF] Manière | Modeling sintering anisotropy in ceramic stereolithography of silica[END_REF].

The "threshold model" considers a critical porosity of 𝜃 𝑐 = 0.068 in the bulk modulus equation [START_REF] Chen | 3D printing of ceramics: A review[END_REF], 𝑃 𝑠 = 0 and: [START_REF] Wang | Review of additive manufacturing methods for high-performance ceramic materials[END_REF].

𝜂 0𝑟 = 𝜂 0𝑧 + 0.00021(𝜃 -0.068) 2.7 (
We can note the porosity expression of 𝜂 0𝑟 also has a critical porosity because the underlying hypothesis of the "threshold model" is the presence of large size and stable pores at the final stage. With such hypothesis the anisotropy should vanish when these stable pores appear.

Simulation of the sintering deflection test

In this section, the "pore gas pressure model" and "threshold model" are tested on the deflection case; the results are shown in figure 4. The two lateral supports correspond to the alumina cup edge represented in figure 1a. In these simulations the stress developed along the bar direction is in tension below the bar and in compression on the upper part. However, the model identified by dilatometry (1X) predicts an excessive deflection compared to the experiment (0.5 mm). This difference is supposed to be attributed to the preferential orientation of the extruded rods along the bar while typical printed samples in figure 2 are a sequence of layers with different rod orientations. In order to slow the creep deformation behavior while keeping the anisotropic sintering behavior unchanged, one way is to modify exclusively the shear creep viscosity by an amplification factor. By doing so, only the shear components are affected without modifying the diagonal term of the tensors component in equation [START_REF] Chartier | Rapid Prototyping of Ceramics[END_REF]. The first observation of figure 4 is that a high amplification factor is required to significantly slow down the shear deformation. In the "pore gas pressure model" the effect is very limited, only the "threshold model" seems to efficiently react to this additional factor. In order to find the optimum shear factor, the modeled deflection results have been plotted in figure 5a. The optimum factor is 44, the corresponding simulation is reported in figure 5b with (as comparison) the lateral view of the bar deflection test showing similar deflection levels. At this stage we obtain a model reproducing the anisotropic sintering response and the high creep resistance of preferential rods orientation in bars shapes. However, the reason for this high creep resistance is questioning and in the next section, we will explore two possible reasons for this creep resistance by mesoscale simulations.

Mesoscale sintering simulation

Before simulating the creep behavior at the microstructure level, a careful analysis of the sintered porcelain microstructure is required. The optical microscopy and EDS mapping of the polished sample are reported in figure 6a andb. The microstructure is homogeneous as shown in figure 6c. Vitrification is the mechanism by which porcelain sintering occurs: the microstructure typically shows an important proportion of glass phase with multimodal porosity size and the presence of large silica particles trapped in the glassy phase. In the EDS map of Al, Si elements on figure 6b, the silica phase appears in green, the mixed violet/green area is the glassy phase, and black zones are porosity. As highlighted in the previous paper [START_REF] Manière | 3D printing of porcelain and finite element simulation of sintering affected by final stage pore gas pressure[END_REF], the sintering of this material suffers from a swelling phenomenon which induces an important pore growth in the interlayer zone. From this microstructure observation, two origins for the creep resistance will be explored. On one hand, the align rods and the close pore gas pressure at the interlayer may implicate a higher effective consistency in shear solicitation. On the other hand, the silica large particles preferentially oriented in the extrusion direction (see figure 6b) may also increase the effective shear viscosity by playing a similar role than rigid inclusions.

Upper

Middle

Lower

These two approaches will be tested in a mesoscale simulation with the aim of testing the magnitude of both hypotheses on the effective viscosity amplification. As represented in figure 7a, the mesoscale simulation is built as a virtual shear test at microstructure level. For neglecting the densification effect, the porosity of the bulk is taken to the value of 0.1, a value close to the overall porosity at the beginning of the holding where the deflection phenomenon is active. In this shear configuration, the shear creep strain rate expression is:

𝜀̇𝑦 𝑧 = 1 2𝜂 𝜎 𝑦𝑧 𝜑 (14)
.

The sliding velocity 𝑣 𝑔 is related to the strain rate by: [START_REF] Du | 3D printing of ceramic-based scaffolds for bone tissue engineering: an overview[END_REF].

𝜀ẏ 𝑧 = |𝑣 𝑔 | 2𝑒 (
The effective viscosity is then calculated by: [START_REF] Miyanaji | Process development for green part printing using binder jetting additive manufacturing[END_REF].

𝜂 = 𝑒 |𝑣 𝑔 | 𝜎 𝑦𝑧 𝜑 (
In the simulation illustrated in figure 7, the effect of the interlayer higher porosity is tested.

The thickness of the layer and interlayer porous zone has been determined by the microstructure on figure 6a. The shear displacement is shown in figure 7b ; the term -P l +P gas can be assimilated to the effective sintering stress. The latter has a negative value when the curvature capillarity forces dominate (densification) and a positive value when trap gas pressure dominates (swelling). With the chosen porosity, the mesoscale simulation (figure 7b) predicts a situation of swelling in the interlayer region and densification in the layer bulk.

This situation is very close to the microstructure seen on figure 6 with large pores at the interfaces. This particular situation allows the layer to slide over each other and generates a highly effective shear viscosity more than ten times the bulk viscosity (figure 7c). The local maximum at 1200 s in figure 7c seems to correspond to the onset of intense swelling (P l <<P gas and significant increase of the porous layer thickness) . In figure 8a, the presence of SiO 2 rigid inclusion is tested. The simulation with 25% of SiO 2 was based on the microstructure illustrated on figure 6a. Higher contents are also tested. The shear effective viscosity is reported in figure 8b, with in figure 8c the shear deformations. In this approach, we assume that the bulk viscosity is similar to the effective viscosity of the glassy phase responsible for sintering densification. For the 25%, the amplification of the effective viscosity is only twice the bulk. However, assuming the SiO 2 is partially dissolved into the glassy phase during sintering, the zone of rigid influence may be higher. Therefore a higher content of SiO 2 may also be considered even if their amplification influence on the effective viscosity is only five times the bulk viscosity. This mesoscale simulation shows the effective viscosity for creep deformation may be higher than the effective viscosity calculated from the sintering phenomenon. Both SiO 2 inclusions and interlayer structure swelling may play a role in the creep resistance experimentally observed. The first aspect to take into account for modeling the cup shapes sintering from green state is the drying step that modifies the 3D design dimensions. Without a 3D scanner or a numerical simulation tool able to predict the drying stage, a correction on the height and radial dimension was applied to match the green shape dimensions from the printed CAD (Computer Aid Design) shapes. The dimensions of the green printed cups were measured accordingly to the scheme on figure 9a (r 1 ,r 2 ,r 3 ,r 4 and h). In figure 9bc, the radial dimensions of the cup were fitted with a polynomial function allowing applying the radial dry shape correction to the FEM green shape. In the same way, the vertical dimensions of the CAD file were corrected based on the cup height information. This approach is not ideal as not all the dry radial dimensions matched (only two of them are approached).

Sintered

The results of the simulation of the small and large cups by "pore gas pressure" and "threshold" approach are illustrated in figure 10ab. For the "threshold" approach, the shear viscosity amplification factor [START_REF] Manière | Porosity dependence of powder compaction constitutive parameters: Determination based on spark plasma sintering tests[END_REF] was taken into account for the lateral parts of the shape (without the support) because this zone has a preferential orientation like for the bar shape (figure 4b). The simulation figure 10ab shows the stress developed in the shape by the action of its own weight. As expected, the stress is higher for the large cup. The internal stress is concentrated in the zone "r2" (see figure 9a). This zone has a cantilever shape and supports the weight of the upper part. If the wall thickness had been thinner in this area (like for lattice structures), a distortion could have happened.

In figure 10c, the experimental/simulated dimensions of the green sintered and simulated sintered shapes are compared. The ratio of the final and initial dimensions (df/di) highlights the anisotropic behavior. The dimensions of the two simulation approaches are very close; the error is in most cases below 5%. For the small cup, a higher error close to 8% is noted and can be explained by the lower experimental anisotropy of this small cup shape. It is possible that the thin wall of this shape implies a different extruded rod distribution and then a different sintering behavior.

Conclusion

In this study, a comprehensive simulation method of robocasted ceramic sintering has been presented. This approach encompasses the identification of the anisotropic sintering behavior, the adjustment of the sintering model via deflection sintering tests, and the finite element simulation of a complex shape. The aim of this study is to establish a methodology able to simulate the sintering of 3D printed ceramic object. This numerical tool is required to predict at the design stage the shape distortions of the sintering stage. This modeling approach is aimed at helping the design of the shapes and supports at the conception stage for production purpose. The following points are highlighted:

 The anisotropic sintering model can be identified by the adjustment of dilatometry and deflection sintering tests. The sintering behavior and creep deformation at high temperature are determined experimentally by this method.

 The model can take into account late stage sintering swelling phenomena and the shrinkage anisotropy originated from a heterogeneous porosity distribution during the printing.

 The printed specimens have a layered structure implying a noticeable creep deformation resistance at high temperature. The deflection test showed a very limited deformation. From mesoscale simulations, this phenomenon may be explained by a combination of SiO 2 particle acting like rigid inclusions and the trap gas pressure in the interlayer area generating a localized layered stress field.

 The simulation tool predicts the sintering shrinkage with precision below 5%. A discrepancy is noted on printed shapes with a thin wall probably due to a different printing layered structure at the edges.

Figure captions

Table captions

Table 2 Simulation/experiment dimensions comparison

  Porosity rate (s -1 ) 𝜎 Stress tensor (N.m -2 ) 𝜎 𝑒𝑞 Equivalent stress (N.m -2 ) 𝜀̇ Strain rate tensor (s -1 ) 𝜀ė 𝑞 Equivalent strain rate (s -1 ) 𝜀ṙ Radial strain rate component (s -1 ) 𝜀ż Z axis strain rate component (s -1 ) 𝜀ṡ 𝑖𝑛𝑡𝑒𝑟𝑖𝑛𝑔 Free strain rate due to sintering (s -1 ) 𝜑 Shear modulus 𝜓 Bulk modulus 𝜃 𝑐 Critical porosity for sintering final stage Pl Sintering stress (Pa) 𝑃 𝑠 Pore gas pressure in closed porosity (Pa) 𝕚 Identity tensor 𝛼 Surface energy (J.m -2 ) 𝑟 Grain radius (m) 𝑒̇ Trace of the strain rate tensor (s -1 ) 𝜂 Material viscosity (Pa.s) 𝜂 𝑟 Radial component of the viscosity (Pa.s) 𝜂 𝑧 Z component of the viscosity (Pa.s) 𝜂 0 Viscosity pre-exponential factor (Pa.s) 𝜂 0𝑟 Radial component of the viscosity pre-exponential factor (Pa.s) 𝜂 0𝑧 Z component of the viscosity pre-exponential factor (Pa.s) 𝑄 Viscosity activation energy (J.mol -1 ) R Gas constant 8.314 (J.mol -1 .K -1 ) T Temperature (K) 𝜎 𝑦𝑧 YZ stress component (N.m -2 )

Figure 1

 1 Figure 1 Bar deflection sintering test, a) configuration before sintering, b) lateral view of the sintered specimen.

Figure 2

 2 Figure 2 Sintering anisotropy parametric analysis: a) effect of heating rate of 10 mm specimens, b) scale effect on specimens with diametric homothetic increase from 10 to 30 mm, c) filling ratio effect on 30 mm diameter samples, the filling architecture is shown on lower right by the yellow walls in the volume of the parts.

Figure 3

 3 Figure 3 Sintering model identification from dilatometry tests, a) adjustment of the modeled specimen diameter and height, b) experimental and modeled porosity.

Figure 4 Figure 5 a

 45 Figure 4 Modeling of the sintering deflection test by pore gas (a) and threshold (b) approaches, the creep based deflection magnitude has been modified via a parametric study on the shear viscosity terms.

Figure 6 a

 6 Figure 6 a) Optical microscopy of the sintering specimen microstructure, b) Al, Si elements EDSmap of the latter microstructure c) microstructure at the upper, middle, and lower part of a 10 mm diameter sintered specimen.

Figure 7

 7 Figure 7 Mesoscale simulation of the layered microstructure assuming pore gas pressure and lower porosity in the layers interfaces, a) microstructural model structure, b) mesoscale simulated images of the shear deformations and effective sintering stress (-P l +P gas ), c) simulated apparent shear viscosity vs continuum viscosity.

Figure 8 4 . 5 .

 845 Figure 8 Mesoscale simulation of the layered microstructure assuming SiO 2 rigid inclusions, a) microstructural model structure, b) simulated apparent shear viscosity vs continuum viscosity for different content of SiO 2 , c) mesoscale simulated images of the shear deformations.

Figure 9 a

 9 Figure 9 a) Cup shape main reference dimensions, b,c) experimental radial corrective factor vs cup height to assumes the shape drying shrinkage on the initial simulated geometry for the small cup (b) and large cup (c).

Figure 10

 10 Figure 10 Simulated small and large cup sintering with pore gas model (black lines correspond to the green shape dimensions) (a) and the threshold approach with shear viscosity correction (b).

Figure 1

 1 Figure 1 Bar deflection sintering test, a) configuration before sintering, b) lateral view of the

Figure 2

 2 Figure 2 Sintering anisotropy parametric analysis: a) effect of heating rate of 10 mm

Figure 3

 3 Figure 3 Sintering model identification from dilatometry tests, a) adjustment of the modeled

Figure 4

 4 Figure 4 Modeling of the sintering deflection test by pore gas (a) and threshold (b)

Figure 5 a

 5 Figure 5 a) Deflection curves for different shear viscosity factors, b) experimental vs optimum

Figure 6 a

 6 Figure 6 a) Optical microscopy of the sintering specimen microstructure, b) Al, Si elements

Figure 7

 7 Figure 7 Mesoscale simulation of the layered microstructure assuming pore gas pressure and

Figure 8

 8 Figure 8 Mesoscale simulation of the layered microstructure assuming SiO 2 rigid inclusions,

Figure 9 a

 9 Figure 9 a) Cup shape main reference dimensions, b,c) experimental radial corrective factor

Figure 10

 10 Figure 10 Simulated small and large cup sintering with pore gas model (black lines

Table 1 Simulation/experiment dimensions comparison Green body Sintered EXP df / di EXP Sintered SIM with pore gas Sintered SIM with SHEAR

 1

		h	32.74	28.24	0.86	25.78	25.92
		r 1	9.88	8.68	0.88	8.74	8.75
	Small Cup	r 2 r 3	18.18 20.66	16.16 18.53	0.89 0.90	16.38 18.10	16.39 18.12
		r 4	23.32	20.82	0.89	20.70	20.71
		h	60.74	49.23	0.81	47.80	47.38
		r 1	18.41	15.90	0.86	16.30	16.24
	Large Cup	r 2 r 3	32.81 37.21	29.40 33.30	0.90 0.90	29,82 32.75	29.7 32.64
		r 4	42.22	38.14	0.90	37.69	37.54

Acknowledgements

The help and support of Jérôme Lecourt and Christelle Bilot is gratefully acknowledged.

Funding

No funding

Data availability

All data is published with the paper.

Compliance with ethical standards

Conflict of interest: The authors declare that they have no conflict of interest.

Consent to participate: No human or animal was involved in this work; thus, no consent was required.

Consent to publish: All authors have given their permission for publishing this work.

Code availability:

The simulation software Comsol Multiphysics 5.5 is a commercially licensed product of Comsol Corporation.


Charles Manière
email: charles.maniere@ensicaen.fr

Christelle Harnois

Guillaume Riquet

Jérôme Lecourt

Christelle Bilot

Sylvain Marinel

Flash spark plasma sintering of zirconia nanoparticles: electro-thermalmechanical-microstructural simulation and scalability solutions

Keywords: Flash sintering, zirconia, simulation, scalability, spark plasma sintering

Flash sintering is a very promising method with a great potential for the ultra-rapid production of objects with firing processes of mere seconds. However, the transition from few millimeters samples to larger scales is a key issue. In this study, we explore the spark plasma sintering approach allowing a stable hybrid heating in flash condition and producing near net shape sintered specimens. Two electric current configurations are employed for diameters of 15 and 30 mm to determine the effect of scale change and of electric current concentration.

These experiments coupled with a Multiphysics simulation indicate that stable flash conditions can be reached for 30 mm specimens. However, even if the electrical current concentration is very effective at small sizes, it generates peripheral hot spots for large specimen dimensions. The blackening effect on zirconia flash specimens acts then as an indicator of the specimen thermal history.

Conventional SPS

Flash SPS

Introduction

Since the discovery of flash sintering in 2010 [START_REF] Cologna | Flash Sintering of Nanograin Zirconia in <5 s at 850°C[END_REF][START_REF] Raj | Methods of flash sintering[END_REF], many efforts have been employed to study the underlying heating/sintering mechanisms and to scale up this ultra-rapid process [START_REF] Biesuz | Flash sintering of ceramics[END_REF][START_REF] Dancer | Flash sintering of ceramic materials[END_REF][START_REF] Yu | Review of flash sintering: materials, mechanisms and modelling[END_REF].

The flash sintering involves an abrupt heat dissipation followed by the sintering of the specimen in mere seconds. For zirconia, this "flash event" is mostly initiated by the NTC (negative temperature coefficient) resistivity of this material generating an incubation stage and an abrupt current flowing in the ceramic at high temperature [START_REF] Todd | Electrical characteristics of flash sintering: thermal runaway of Joule heating[END_REF]. The ultra-rapid sintering kinetics are generally not explained by the heating speed that would require higher temperatures and suggest additional mass transport mechanism or kinetic mechanisms [START_REF] Cologna | Flash Sintering of Nanograin Zirconia in <5 s at 850°C[END_REF][START_REF] Raj | Influence of Externally Imposed and Internally Generated Electrical Fields on Grain Growth, Diffusional Creep, Sintering and Related Phenomena in Ceramics[END_REF].

Concerning the non-thermal mass transport mechanisms, we can cite the current effects (electromigration [START_REF] Trzaska | Electromigration experiments by spark plasma sintering in the silver-zinc system[END_REF][START_REF] Olevsky | Constitutive modeling of spark-plasma sintering of conductive materials[END_REF], electroplasticity [START_REF] Conrad | Electroplasticity in metals and ceramics[END_REF][START_REF] Lee | Electric current effects in spark plasma sintering: From the evidence of physical phenomenon to constitutive equation formulation[END_REF][START_REF] Lee | Effect of electric current on densification behavior of conductive ceramic powders consolidated by spark plasma sintering[END_REF], Frenkel pair cascade [START_REF] Raj | Influence of Externally Imposed and Internally Generated Electrical Fields on Grain Growth, Diffusional Creep, Sintering and Related Phenomena in Ceramics[END_REF]) and the field effects (for instance, the ponderomative forces [START_REF] Rybakov | The microwave ponderomotive effect on ceramic sintering[END_REF][START_REF] Olevsky | Ponderomotive effects during contact formation in microwave sintering[END_REF][START_REF] Chaim | Electric field effects during spark plasma sintering of ceramic nanoparticles[END_REF]). Concerning the thermal/kinetics effects, we can cite the localized heating [START_REF] Raj | Influence of Externally Imposed and Internally Generated Electrical Fields on Grain Growth, Diffusional Creep, Sintering and Related Phenomena in Ceramics[END_REF][START_REF] Narayan | A new mechanism for field-assisted processing and flash sintering of materials[END_REF][START_REF] Chaim | On thermal runaway and local endothermic/exothermic reactions during flash sintering of ceramic nanoparticles[END_REF], the high heating rate effect on delayed pore surface diffusion [START_REF] Olevsky | Consolidation enhancement in spark-plasma sintering: Impact of high heating rates[END_REF][START_REF] Ji | Ultra-fast firing: Effect of heating rate on sintering of 3YSZ, with and without an electric field[END_REF]. We may also cite the non-equilibrium state of grain boundaries [START_REF] Zhang | The nature of grain boundaries in alumina fabricated by fast sintering[END_REF][START_REF] Chaim | Particle Surface Softening as Universal Behaviour during Flash Sintering of Oxide Nano-Powders[END_REF] and grain boundary complexion [START_REF] Cantwell | Grain boundary complexions[END_REF][START_REF] Cantwell | Expanding timetemperature-transformation (TTT) diagrams to interfaces: A new approach for grain boundary engineering[END_REF]. The origin of the flash sintering ultra-rapid mechanism is still under hard debate [START_REF] Bordia | Current understanding and future research directions at the onset of the next century of sintering science and technology[END_REF][START_REF] Olevsky | Field-Assisted Sintering[END_REF]. It must bear in mind that the high thermal instability of flash sintering makes the estimation of the effective temperature very difficult due to the presence of hot spots and temperatures gradients [START_REF] Biesuz | Flash sintering of ceramics[END_REF][START_REF] Jones | Promoting microstructural homogeneity during flash sintering of ceramics through thermal management[END_REF][START_REF] Qin | Temperature gradient and microstructure evolution in AC flash sintering of 3 mol% yttria-stabilized zirconia[END_REF]. It is then difficult to separate the thermal and non-thermal effects.

In flash sintering, the cooling flux generates thermal gradient at the edges of the specimen [START_REF] Manière | Impact of convection and radiation on direct/hybrid heating stability of field assisted sintering[END_REF]. The NTC behavior tends to amplify these thermal gradients by a higher volumetric dissipation in the hot area where the current tends to concentrate [START_REF] Manière | Inherent heating instability of direct microwave sintering process: Sample analysis for porous 3Y-ZrO2[END_REF][START_REF] Manière | Fully coupled electromagnetic-thermalmechanical comparative simulation of direct vs hybrid microwave sintering of 3Y-ZrO2[END_REF]. For materials having low thermal conductivity like zirconia, this phenomenon is even more amplified and hot spots may appear. Dong [START_REF] Dong | On the Hotspot Problem in Flash Sintering[END_REF] has expressed this in a figure of merit based on thermal relaxation theory. In microwave sintering, the hot spot formation is typically observed in the case of the direct heating of zirconia [START_REF] Manière | Inherent heating instability of direct microwave sintering process: Sample analysis for porous 3Y-ZrO2[END_REF]. The latter is initiated by cooling fluxes by conduction, convection and radiation at the edge [START_REF] Manière | Impact of convection and radiation on direct/hybrid heating stability of field assisted sintering[END_REF]. The finite element simulation helps in predicting and quantifying these fluxes. For SPS [START_REF] Risbud | Preface and historical perspective on spark plasma sintering[END_REF], the cooling fluxes at the specimen boundaries are only caused by the conduction through the surrounding graphite elements acting as susceptors in direct contact. The thermal gradients are strongly influenced by the electrical and thermal contacts which play a similar role are thermal barrier and dissipate heat at the boundaries [START_REF] Manière | Contact resistances in spark plasma sintering: From in-situ and ex-situ determinations to an extended model for the scale up of the process[END_REF][START_REF] Manière | Finiteelement modeling of the electro-thermal contacts in the spark plasma sintering process[END_REF][START_REF] Wei | Experimental Investigation of Electric Contact Resistance in Spark Plasma Sintering Tooling Setup[END_REF][START_REF] Zavaliangos | Temperature evolution during field activated sintering[END_REF]. These hybrid heating characteristics concentrate the heat in the powder making flash SPS more stable and more controllable than the conventional resistive approach. In addition, the punches and die allow controlling the external shape of the flashed specimen.

Different flash SPS approaches have been developed. The first is based on a sinter-forging configuration where a graphite felt or foil is used to initiate the heating. Zirconia [START_REF] Vasylkiv | Flash spark plasma sintering of ultrafine yttria-stabilized zirconia ceramics[END_REF],

ZrB 2 [START_REF] Grasso | Flash Spark Plasma Sintering (FSPS) of Pure ZrB 2[END_REF], SiC [START_REF] Grasso | Flash Spark Plasma Sintering (FSPS) of α and β SiC[END_REF][START_REF] Olevsky | Flash (Ultra-Rapid) Spark-Plasma Sintering of Silicon Carbide[END_REF][START_REF] Grasso | Ultra-Rapid Crystal Growth of Textured SiC Using Flash Spark Plasma Sintering Route[END_REF][START_REF] Zapata-Solvas | Preliminary investigation of flash sintering of SiC[END_REF], Nd-Fe-B [START_REF] Castle | Rapid sintering of anisotropic, nanograined Nd-Fe-B by flash-spark plasma sintering[END_REF], TiB 2 -hBN [START_REF] Mckinnon | Flash spark plasma sintering of cold-Pressed TiB 2 -h BN[END_REF], B 4 C [START_REF] Niu | Ultra-fast densification of boron carbide by flash spark plasma sintering[END_REF] materials have been produced by this method. This method is very interesting for achieving high density of refractories as the radial shrinkage allows the pore flattening and porosity elimination during the sinterforging SPS approach [START_REF] Tünçay | Spark plasma sintering and spark plasma upsetting of an Al-Zn-Mg-Cu alloy[END_REF][START_REF] Shang | Evolution of anisotropy in hierarchical porous ceramics during sinter-forging[END_REF]. High dimension specimens have been produced by this method [START_REF] Grasso | Flash Spark Plasma Sintering (FSPS) of α and β SiC[END_REF]. Another approach uses a graphite die to control the final shape of the sample.

The graphite die is generally electrically insulated by a boron nitride spray to concentrate the electrical current in the powder and the graphite foil which initiate the heating and act as a susceptor. This configuration has been used to sinter ZrB 2 , MoSi 2 , Al 2 O 3 [START_REF] Zapata-Solvas | Ultra-fast and energy-efficient sintering of ceramics by electric current concentration[END_REF], SiC [START_REF] Román-Manso | Effects of Current Confinement on the Spark Plasma Sintering of Silicon Carbide Ceramics[END_REF]. In a previous study [START_REF] Manière | All-Materials-Inclusive Flash Spark Plasma Sintering[END_REF], it was shown that this configuration provides an electro-thermalmechanical confinement to the powdered specimen that can be ultra-rapidly heated and sintered with very controllable and homogeneous conditions. In the latter, we showed that these ultra-rapid conditions allowed the sintering of almost all materials from dielectrics, semiconductors to metals. These experiments were done on alumina, zirconia and nickel and for 10 mm diameter samples. However, the preservation of these stable ultra-rapid conditions for bigger samples is still questionable.

In this work, the scalability of flash SPS of zirconia is tested for cylindrical samples with diameter ranging from 15 to 30 mm. The flash sintering is investigated in different configurations of electrical current to identify stable heating/sintering conditions. A fully coupled electro-thermal-mechanical-microstructural simulation of the different tests is conducted to model the temperature, porosity and grain size fields of each flash test.

Theory and calculations

Electro-Thermal-Mechanical-Microstructural model (ETMM Model)

The ETMM simulation of spark plasma sintering includes the coupled Joule heating and sintering governing equations. The Joule heating obeys the charge conservation and heat transfer equations.

∇. 𝐽 ⃗ = ∇. (𝜎 𝑒𝑙𝑒𝑐 𝐸 ⃗⃗ ) = 0 (1)
∇. (-𝜅∇𝑇) + 𝜌𝐶 𝑝 𝑑𝑇 𝑑𝑡 = 𝑱𝑬 (2)
Where 𝑱𝑬 is the volumetric electrical power dissipated.

The sintering part is simulated with the continuum theory of sintering [START_REF] Bordia | Current understanding and future research directions at the onset of the next century of sintering science and technology[END_REF][START_REF] Olevsky | Theory of sintering: from discrete to continuum[END_REF][START_REF] Olevsky | Fundamental Aspects of Spark Plasma Sintering: II. Finite Element Analysis of Scalability[END_REF]. The local momentum equation of the mechanical part is:

∇. 𝜎 = 0 (3).
The porous viscoplastic sintering behavior is defined by the strain rate and stress tensor relationship:

𝜎 = 𝜎 𝑒𝑞 𝜀̇𝑒 𝑞 (𝜑𝜀̇+ (𝜓 - 1 3 𝜑) 𝑡𝑟(𝜀)𝕚) + 𝑃𝑙𝕚 (4).
The expression of the shear 𝜑 and bulk 𝜓 moduli (porosity functions) can be theoretically approximated [START_REF] Olevsky | Theory of sintering: from discrete to continuum[END_REF][START_REF] Skorohod | Rheological basis of the theory of sintering[END_REF]. However, a SPS-based experimental identification method has been developed [START_REF] Manière | Identification of the Norton-Green Compaction Model for the Prediction of the Ti-6Al-4V Densification During the Spark Plasma Sintering Process[END_REF][START_REF] Martins | In-situ creep law determination for modeling Spark Plasma Sintering of TiAl 48-2-2 powder[END_REF][START_REF] Manière | Porosity dependence of powder compaction constitutive parameters: Determination based on spark plasma sintering tests[END_REF]. For the ceramics the independent identification of the moduli and creep parameters is very challenging due to the grain growth disturbances on the creep behavior at the final stage. A method allowing the direct identification of both the creep and moduli directly from porous stage with tests like sinter-forging has been developed [START_REF] Manière | Porous stage assessment of pressure assisted sintering modeling parameters: a ceramic identification method insensitive to final stage grain growth disturbance[END_REF]. Below, moduli and creep law was obtained for the zirconia powder employed in this study:

𝜑 = (1 - 𝜃 0.6 ) 3.4
(5)

𝜓 = 0.07 (1-𝜃) 5.8 𝜃 (6)
𝜀̇𝑒 𝑞 = 1.34𝐸3[𝑃𝑎 -2 𝑠 -1 𝐾] 𝑇 ( 𝐺 0 𝐺 ) 𝑚 exp ( -536[𝑘𝐽 𝑚𝑜𝑙 -1 ] 𝑅𝑇 )𝜎 𝑒𝑞 2 (7)
with the expressions:

𝜀ė 𝑞 = 1 √1-𝜃 √𝜑𝛾̇2 + 𝜓𝑒̇2 (8)
and the strain rate tensor invariants:

𝑒̇= 𝜀ẋ + 𝜀ẏ+𝜀ż (9) 𝛾̇= √2(𝜀̇𝑥 𝑦 2 + 𝜀ẋ 𝑧 2 +𝜀ẏ 𝑧 2 ) + 2 3 (𝜀ẋ 2 +𝜀ẏ 2 +𝜀ż 2 ) - 2 3
(𝜀ẋ𝜀ẏ + 𝜀ẋ𝜀ż + 𝜀ẏ𝜀ż) [START_REF] Conrad | Electroplasticity in metals and ceramics[END_REF].

The volume change rate is related to the porosity elimination rate by the mass conservation equation: [START_REF] Lee | Electric current effects in spark plasma sintering: From the evidence of physical phenomenon to constitutive equation formulation[END_REF].

𝜃 (1-𝜃) = 𝜀̇𝑥 + 𝜀̇𝑦+𝜀̇𝑧 (
The grain growth is modeled by the following rate equation.

Ġ = 𝐾(𝑇) 𝐺 𝑝 (12)
𝐾(𝑇) = 𝐾 0 exp ( -𝑄 𝐺 𝑅𝑇 ) (13)

Boundary and interface conditions

The electro-thermal boundary conditions and mesh of the SPS column are represented in figure 1. An anisotropic mesh is employed to have three elements in the 0.2 mm thickness of the graphite foil where the electrical current is concentrated in flash conditions. The electrode water cooling is modeled by a convective flux [START_REF] Manière | Finiteelement modeling of the electro-thermal contacts in the spark plasma sintering process[END_REF][START_REF] Pavia | Electro-thermal measurements and finite element method simulations of a spark plasma sintering device[END_REF][START_REF] Molénat | Temperature Control in Spark Plasma Sintering: An FEM Approach[END_REF] with the following equation:

𝜑 𝑐𝑜𝑛𝑣 = ℎ 𝑐𝑜𝑛𝑣 (𝑇 𝑤𝑎𝑡𝑒𝑟 -𝑇 𝑠𝑜𝑙𝑖𝑑 ) (14)
.

The graphite surfaces radiate towards each other and metallic adjacent surfaces of the SPS are maintained at room temperature. The heat loss by radiation is calculated by the surface to surface radiation module of comsol defining the radiosity flux 𝐽 𝑟𝑎𝑑 as the sum of the surface thermal radiation ϵ𝑒 𝑏 (𝑇) and the reflected part (refl) of the incident flux G.

𝐽 𝑟𝑎𝑑 = 𝑟𝑒𝑓𝑙 + ϵ𝑒 𝑏 (𝑇) = (1 -ϵ)𝑮 + ϵ(𝑁𝑟) 2 𝜎 𝑠 𝑇 4 (15)
.

For a surface, the net inward heat flux 𝜑 𝑟𝑠𝑠 is defined by:

𝜑 𝑟𝑠𝑠 = ϵ(𝑮 -𝑒 𝑏 (𝑇)) (16)
.

Figure 1 Flash SPS configuration, (a) 30 mm samples configuration and model boundary condition (b) mesh.

A graphite emissivity of 0.8 is chosen [START_REF] Manière | Contact resistances in spark plasma sintering: From in-situ and ex-situ determinations to an extended model for the scale up of the process[END_REF]. The natural convection [START_REF] Manière | Fluid dynamics thermo-mechanical simulation of sintering: Uniformity of temperature and density distributions[END_REF][START_REF] Manière | Graphite creep negation during flash spark plasma sintering under temperatures close to 2000[END_REF] is neglected since the cavity is under vacuum. The electric boundary conditions are the following: the ground is located at the lower electrode, a numerical PID is incorporated [START_REF] Manière | Proportional integral derivative, modeling and ways of stabilization for the spark plasma sintering process[END_REF]

The electric and thermal contact resistances (ECR and TCR) have a high influence on the thermal field [START_REF] Manière | Finiteelement modeling of the electro-thermal contacts in the spark plasma sintering process[END_REF][START_REF] Wei | Experimental Investigation of Electric Contact Resistance in Spark Plasma Sintering Tooling Setup[END_REF][START_REF] Anselmi-Tamburini | Fundamental investigations on the spark plasma sintering/synthesis process[END_REF][START_REF] Maniere | Pulse analysis and electric contact measurements in spark plasma sintering[END_REF][START_REF] Vanmeensel | Modelling of the temperature distribution during field assisted sintering[END_REF] and obey:

𝐽 𝑐 = 1 𝐸𝐶𝑅 (𝑈 1 -𝑈 2 ) (18) 𝑞ċ = 1 𝑇𝐶𝑅 (𝑇 1 -𝑇 2 ) (19)
.

The ECR TCR and temperature-dependent electro-thermal properties have been taken from [START_REF] Manière | Contact resistances in spark plasma sintering: From in-situ and ex-situ determinations to an extended model for the scale up of the process[END_REF] and [START_REF] Van Der Laan | Simple method for the identification of electrical and thermal contact resistances in spark plasma sintering[END_REF]. The electrical conductivity of zirconia has been taken from [START_REF] Xu | DC electric field-assisted hot pressing of zirconia: methodology, phenomenology, and sintering mechanism[END_REF] where flash sintering was applied with an alumina die without graphite foil, allowing to get an estimation of the zirconia resistivity under vacuum conditions at high temperatures. In the calculation, the boron nitride spray interfaces are assumed to be electrically insulating. All the material properties are gathered in the appendix.

Experiment and method

All experiments have been conducted on the SPS apparatus "FCT HPD25" with the graphite grade mersen 2333. A 40 nm yttria stabilized zirconia nano-powder (Tosoh TZ-3Y-E) is employed for this study. Two tooling geometries are employed. The geometry of the 30 mm sample diameter tests is reported in figure 1a. For the 15 mm tests, the heights of the punches and die are similar and the die diameter is reduced to 45 mm. The sintering sample height is 2 mm for all experiments, 2 and 8 g of powder are taken for the 15 and 30 mm tests, respectively. Two electrical current configurations are used: one conventional (called "die heating") where the current can flow through the punches, graphite foil, die and another referred as "concentrated" where the die is electrically insulated by a boron nitride spray (figure 1a) to concentrate the current in the powder and the graphite foil acting as a susceptor (called "concentrated" configuration). In order to complete this explorative design, two thermal cycles were imposed to the punch temperature, one typical with a 100 K/min ramp and another at 1000 K/min to impose flash but controlled conditions. A dwell is programmed at 1400 °C but the test is manually stopped when the displacement indicates the end of the sintering process. The powder quantity is calculated for 2 mm thick sintered samples. The experimental design includes 8 tests with the combinations of the three following parameters:

of the specimen diameter (Ø15, Ø30, scalability), heating rate (100, 1000 K/min) and electric current configurations ("die heating", "concentrated"). The samples relative density was measured by Archimedes's method and the average grain size was determined on polished surfaces using the linear intercept method with Mendelson's stereological factor (1.56) [START_REF] Mendelson | Average Grain Size in Polycrystalline Ceramics[END_REF].

An ETMM simulation of each of these tests has been conducted to describe the differences in the Multiphysics fields (temperature, porosity, grain size, etc.). The calculations have been done via the simulation software COMSOL Multiphysics® with the heat transfer module for the surface to surface radiation and electro-thermal contact resistance and the nonlinear material module to implement manually the sintering physics. The obtained microstructures have been analyzed by SEM on polished cross-section surfaces.

Results and discussions

ETMM simulation calibration

The simulation study started by conducting a calibration of the thermal field (with in particular the lateral contact resistance which differs from [START_REF] Manière | Contact resistances in spark plasma sintering: From in-situ and ex-situ determinations to an extended model for the scale up of the process[END_REF][START_REF] Van Der Laan | Simple method for the identification of electrical and thermal contact resistances in spark plasma sintering[END_REF]). Once the ECR, TCR have been adjusted with the punch/die pyrometers, the grain growth and the creep grain size exponent m were adjusted to correspond to the final relative density and grain size for all tests. As said earlier, the moduli and creep behavior was determined previously [START_REF] Manière | Porous stage assessment of pressure assisted sintering modeling parameters: a ceramic identification method insensitive to final stage grain growth disturbance[END_REF]. The latter employed electrically insulated and isothermal conditions at the porous stage. The following values were obtained k 0 = 0.04 m 1+p .s -1 , Q G = 575 kJ/mol, p = 1.63 (from [START_REF] Manière | Modeling zirconia sintering trajectory for obtaining translucent submicronic ceramics for dental implant applications[END_REF]) and m = 2. The graphite/foil/graphite ECR and TCR formulas obtained are 3.11E -7 exp(-2.61E - shown in figure 2 where the ETMM simulation was conducted and compared to the experimental values for the test Ø30 mm, 100 K/min, in "die heating" current configuration.

SPS tests curves

The punch temperature, SPS displacement, RMS current (I rms ) and voltage (U rms ) curves are reported in figure 3 for the 30 mm diameter tests. These curves allow comparing the heating/sintering response of the four configurations for current distribution and heating rate.

The curves in figure 3a show the typical response of a SPS test at 100 K/min with linear I rms , U rms curves and apparent sintering temperatures taking place between 1000 and 1400 °C. In the same electric current configuration at 1000 K/min, abrupt current/voltage curves are observed with an overshot of the temperature and a quick stabilization to similar apparent sintering temperatures. If the "concentrated" configuration of electrical current is employed (BN spray), the test at 100 K/min (figure 3c) shows a non-linear current profile compared to the equivalent die heating test (figure 3a). At 1000 K/min (figure 3d) an abrupt current profile without stabilization is observed while the sintering duration is 30 sec between 800 and 1400°C of apparent temperatures. This behavior of the temperature, current/voltage experimental curves is similar to the 15 mm tests.

In order to compare the behavior of the apparent sintering response to all tests, the displacement curves are plotted vs the apparent temperature in figure 4a. When using the concentrated configuration, the sintering zones take place about 150 K lower than it does with the die heating method. For the 15 mm diameter tests, the curves at 1000 K/min are systematically shifted toward higher temperatures. For the 30 mm diameter tests, this shift is smaller for the "die heating" configuration and inverted for the "concentrated" configuration where slightly lower sintering temperature is observed for the highest heating rate. These observations suggest a significant temperature difference between the apparent punch temperature and the effective powder temperature. In the next section, the ETMM simulations of all the 8 tests are fully disclosed to describe the internal fields and their differences. 1000 K/min for sintering duration of about 250 sec, and 30 sec, respectively. These curves also show the impact of the grain growth on the slowing down of the densification curves due to equation [START_REF] Raj | Influence of Externally Imposed and Internally Generated Electrical Fields on Grain Growth, Diffusional Creep, Sintering and Related Phenomena in Ceramics[END_REF].

The very interesting aspect of these curves (figure 4b) is that ultra-rapid sintering rates could be obtained without modifying the sintering model parameters identified at the same applied pressure and in isothermal and electrically insulated conditions in ref [START_REF] Manière | Porous stage assessment of pressure assisted sintering modeling parameters: a ceramic identification method insensitive to final stage grain growth disturbance[END_REF]. To carefully describe this, the simulated temperatures disclosed in the next section must be discussed as well. This suggests experimentally that the measured moduli at porous stage give a better description of the true porous skeleton (with pore shape different than spherical pores). This suggests a high sintering reactivity of powders under flash conditions.

Figure 4 (a) Displacement and (b) simulated relative density curves, "dh" and "conc" refer to the "die heating" and "concentrated" configurations respectively. The 15 mm diameter tests in "die heating" electrical current configuration are reported in figure 5. The heating is concentrated in the punches where a temperature of 1900 °C is observed. In the sample, the temperature is distributed from 1600 °C in the center to 1540 °C at the edge where lower grain size and densification are observed. At 1000 K/min, the temperatures/grain sizes are very similar and despite the high heating rate, the densification is close to the one observed on the 100 K/min test. The sintering seems to follow the high heating rates. An overshoot with temperatures 150 K higher is observed for the test at 1000 K/min. The 15 mm diameter tests in "concentrated" electrical current configuration are reported in figure 6. The temperature distribution is very different from the "die heating" configuration.

Is this case, the electrical current is highly concentrated in the powder/graphite and foil/die interface at low temperatures where the zirconia powder behaves practically like a dielectric.

At high temperatures, the zirconia powder becomes conductive and the current lines are distributed between the powder and graphite foil. A significant heat is then dissipated and confined in the powder zone. The specimen temperature has a vertical gradient with cooler temperatures near the contact zone with the punches, and higher temperatures, grain sizes and relative density near the center. Because of this significant heat concentration in the sample, the apparent temperature measured in the punches is 200 K lower than the sample. At 1000 K/min, the latter is increased to 500 K at the end. The temperature distribution is relatively close to the 100 K/min test but the thermal confinement is bigger due to the significant reduction of time for the heat to diffuse from the specimen zone to nearby graphite tooling. The 30 mm diameter tests in "concentrated" electrical current configuration are reported in figure 8. The temperature field distribution from "die heating" to "concentrated" configuration in 30 mm diameter tests is very similar to the test at 15 mm. A significant heat concentration is still observed despite the flatter specimens shape at 30 mm. The axial thermal gradient behavior of the temperature/relative density/grain size fields is still observed. The differences between the apparent punches temperature and the specimen is of the same order as for the 15 mm sample diameter: 200 K and 500 K at 100 and 1000 K/min respectively. In the next section, these results are compared to the experimental microstructures.

Microstructure analysis

The microstructures of all the 15 and 30 mm diameter tests are reported in figure 9. To investigate the microstructure homogeneity and the effects of the Multiphysics field gradients, the microstructure was analyzed in the center and at the edge of the specimens. We can see that the microstructure homogeneity is very satisfactory for all tests in terms of porosity and grain size. The only exception concerns the 30 mm, 1000 K/min, "concentrated" test where an important grain size and porosity gradient are observed near the edge. In the corresponding simulation in figure 8, the radial grain size gradient is observed but the relative density distribution should be axially distributed. In the photograph of this test a complex distribution of blackened zone is observed in the zirconia cross-section. This blackening phenomenon follows the simulated specimen current density distribution reported below in figure 9, it points out a hot spot near the edge. In particular, a thin black zone is present at mid height. A similar distribution has been observed by simulation the Peltier effect in SPS of colusite sulfur thermoelectric specimens [START_REF] Guélou | Issues and opportunities from Peltier effect in functionally-graded colusites: From SPS temperature modeling to enhanced thermoelectric performances[END_REF]. The latter generates local temperature maximums by the coupling of the thermoelectric flux [START_REF] Van Der Laan | Simple method for the identification of electrical and thermal contact resistances in spark plasma sintering[END_REF] and the thermal contact resistance. However, without the information of the Seebeck coefficient at such high temperatures for zirconia, it is difficult to conclude on the origin of this phenomenon. A competition between the reduction front and the pollution species diffusion may also be present.

Concerning the blackening effect on zirconia samples, the aspect of each sintered sample was described in figure 9. Comparing these aspects to the simulated temperatures, this effect seems to appear for the high temperatures > 1500 °C and preferentially during the "concentrated" configurations where the concentration of the electrical current in the powder seems to accelerate the reduction phenomena. The blackening effect is well known is the zirconia flash sintering, it is associated with a reduction phenomenon appearing near the electrodes and propagating through the other [START_REF] Biesuz | Investigation of Electrochemical, Optical and Thermal Effects during Flash Sintering of 8YSZ[END_REF][START_REF] Janek | Electrochemical blackening of yttria-stabilized zirconia â€“ morphological instability of the moving reaction front[END_REF]. The latter is often accompanied by a significant rise of the electrical conductivity switching from ionic conduction to electronic conduction. For SPS imposing reductive conditions, this effect may be accelerated.

Conclusion

In this study, the flash spark plasma sintering of zirconia is investigated from medium to large size samples. The exploration consisted of comparing the typical SPS configuration in electric current distribution and heating rate (100 K/min) to ultra-rapid heating rates (1000 K/min) and concentrated electrical current configuration. The experiments and electro-thermalmechanical-microstructural simulations of each case show that the thermal distribution of the 15 mm and 30 mm sizes are radial for a classic electrical current and becomes axial when the die is electrically insulated. Submicronic microstructures have been obtained with a good stability between the center and the edge. The exception is the 30 mm test in "concentrated" configuration showing intense grain growth and blackening effect at the edges. This shows the scalability limit of the "concentrated" configuration where the die is electrically insulated.

Processing ultra-rapid heating over 15 mm specimens requires distributing the current in the die to avoid excessive overheating in the graphite foil during the preheating.

The thermal-sintering model shows that the sintering follows the high heating rates to sintering time down to 30 sec when 1000 K/min is employed. This suggests the need to use experimentally determined moduli which would help in predicting the high reactivity of the microstructures during ultra-rapid sintering.

The help and support of Christophe Couder, Jérôme Lecourt and Christelle Bilot is gratefully acknowledged. Figure 9 Polished SEM microstructures of the specimens at 100 K/min, 1000 K/min and for the "die heating" and "concentrated" electrical current configurations.

Figure captions

𝐼 1 G

 1 Porosity elimination rate (s -1 ) 𝜎 Stress tensor (N.m -2 ) 𝜎 𝑒𝑞 Equivalent stress (N.m -2 ) 𝜀̇ Strain rate tensor (s -1 ) 𝜀ė 𝑞 Equivalent strain rate (s -1 ) n Creep law stress exponent A Creep law deformability term (s -1 Pa -n ) 𝐴 0 Creep pre-exponential factor (Ks -1 Pa -n ) m Creep law grain size sensibility exponent 𝑄 Sintering activation energy (J.mol -1 ) R Gas constant 8.314 (J.mol -1 .K -1 ) T Temperature (K) 𝜑 Shear modulus 𝜓 Bulk modulus Pl Sintering stress (Pa) 𝕚 Identity tensor 𝑒̇ Trace of the strain rate tensor (s -1 ) 𝛾̇ Shear strain rate invariant (s -1 ) Trace of the stress tensor (N.m -2 ) 𝐺 ̇ Grain growth rate (m.s -1 ) 𝐺 Grain size (m) 𝐺 0 Initial grain size (m) 𝑝 Grain growth equation grain size exponent 𝐾 Grain growth factor (m 1+p .s -1 ) 𝑘 0 Grain growth pre-exponential factor (m 1+p .s -1 ) 𝑄 𝐺 Grain growth activation energy (J.mol -1 ) SPS Spark Plasma Sintering ETMM Electro-Thermal-Mechanical-Microstructural 𝐽 ⃗ Electric current density (A m -2 ) 𝐸 ⃗⃗ Electric field (V m -1 ) J Electric current density norm (A m -2 ) E Electric field norm (V m -1 ) 𝜎 𝑒𝑙𝑒𝑐 Electric conductivity (S m -1 ) 𝜌 Density (kg m -3 ) 𝐶 𝑝 Heat capacity (J kg -1 K -1 ) T Temperature (K) 𝜅 Thermal conductivity (W m -1 K -1 ) 𝜑 𝑐𝑜𝑛𝑣 Convective heat flux (W m -2 ) ℎ 𝑐𝑜𝑛𝑣 Convective heat flux coefficient (W m -2 K -1 ) 𝑇 𝑤𝑎𝑡𝑒𝑟 Water temperature (300 K) 𝑇 𝑠𝑜𝑙𝑖𝑑 Calculated solid temperature (K) 𝐽 𝑟𝑎𝑑 Surface radiosity (W m -2 ) 𝜎 𝑠 Stefan Boltzmann constant (5.67E-8 W m -2 K -4 ) ϵ Emissivity Thermal irradiation flux (W m -2 ) Nr Refractive index 𝑒 𝑏 (𝑇) Surface radiation produced (W m -2 ) 𝑟𝑒𝑓𝑙 Reflected radiative heat flux (W m -2 ) 𝜑 𝑟𝑠𝑠 Net inward radiative heat flux (W m -2 ) 𝑛 ⃗⃗ Normal unit vector 𝐽 𝑛 Imposed normal electric current density (A m -2 ) 𝐽 𝑐 Contact current density (A m -2 ) ECR Electric surface contact resistance (Ω m 2 ) U i Contact i face electric potential (V) 𝑞̇𝑐 Contact heat flux (W m -2 ) TCR Thermal surface contact resistance (m 2 K W -1 ) T i Contact i face temperature (K) α Coefficient of thermal expansion (T -1 )

  at the upper electrode via a normal current density 𝐽 𝑛 and electrical insulations are considered for all other surfaces 𝐽 𝑛 = 0 : . 𝐽 ⃗ = 𝐽 𝑛

Figure 2

 2 Figure 2 Calibration/verification of the electro-thermal-mechanical-microstructural simulation: (a) temperature curves, (b) temperature field of the SPS column, (c) grain size, (d) displacement, (e) specimen final relative density field, (f) simulated average relative density curve.

Figure 3

 3 Figure 3 Temperature, displacement, rms current and voltage experimental data for the 30 mm tests, (a) die heating configuration 100 K/min, (b) die heating configuration 1000 K/min, (c) electric current concentration configuration (concentrated configuration) at 100 K/min, (d) electric current concentration configuration at 1000 K/min.

3 .

 3 Electro-thermal-mechanical-microstructural simulation of the concentrated / die heating configurationsIn this section, the ETMM simulations of all the tests are displayed and compared to each other. Each figure compares the experimental/simulated temperature curves in the punch/die and in the sample, the SPS column simulated temperature field at the end of the heating (the minimum temperature scale is fixed at 900 °C) and the temperature, grain size, relative density of the sample at the end of the cycle.

Figure 5

 5 Figure 5 Simulated ETMM model of the 15 mm diameter specimens in die heating configuration at 100 K/min and 1000 K/min.

Figure 6 Figure 7

 67 Figure 6 Simulated ETMM model of the 15 mm diameter specimens in "concentrated" configuration at 100 K/min and 1000 K/min. The 30 mm diameter tests in "die heating" electrical current configuration are reported in

Figure 8

 8 Figure 8 Simulated ETMM model of the 30 mm diameter specimens in "concentrated" configuration at 100 K/min and 1000 K/min.

Figure 9

 9 Figure 9 Polished SEM microstructures of the specimens at 100 K/min, 1000 K/min and for the "die heating" and "concentrated" electrical current configurations.

Figure 1 Flash

 1 Figure captionsFigure 1 Flash SPS configuration, (a) 30 mm samples configuration and model boundary

Figure 2

 2 Figure 2 Calibration/verification of the electro-thermal-mechanical-microstructural

Figure 3

 3 Figure 3 Temperature, displacement, rms current and voltage experimental data for the

Figure 4 (

 4 Figure 4 (a) Displacement and (b) simulated relative density curves, "dh" and "conc" refer to

Figure 5

 5 Figure 5 Simulated ETMM model of the 15 mm diameter specimens in die heating

Figure 6

 6 Figure 6 Simulated ETMM model of the 15 mm diameter specimens in "concentrated"

Figure 7

 7 Figure 7 Simulated ETMM model of the 30 mm diameter specimens in die heating

Figure 8

 8 Figure 8 Simulated ETMM model of the 30 mm diameter specimens in "concentrated"

T) Ωm 2 and 1.79E -3 exp(-1.53E -3 T) m 2 KW -1 respectively. The outcome of this calibration is

Acknowledgements

This study was funded by the French National Research Agency (ANR), project ULTRARAPIDE N°ANR-19-CE08-0033-01.

Data availability

The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations.

Credit authorship contribution statement

Charles Manière: Conceptualization, Supervision, Modeling, Writing; Christelle Harnois: Conceptualization, Supervision, review & editing; Guillaume Riquet: Review & editing, Christelle Bilot: Characterizations, Jérôme Lecourt: Experimentation, Sylvain Marinel: Conceptualization, Supervision, review & editing.

Appendix

Table A Temperature dependent material properties of graphite, electrodes, and zirconia [START_REF] Manière | Contact resistances in spark plasma sintering: From in-situ and ex-situ determinations to an extended model for the scale up of the process[END_REF][START_REF] Manière | Proportional integral derivative, modeling and ways of stabilization for the spark plasma sintering process[END_REF][START_REF] Manière | Simultaneous spark plasma sintering of multiple complex shapes[END_REF][START_REF] Hayashi | Thermal expansion coefficient of yttria stabilized zirconia for various yttria contents[END_REF].

Materials Expression

Graphite Cp (J .kg -