A Abaza

J Laurencin
email: jerome.laurencin@cea.fr

A Nakajo

S Meille

J Debayle

D Leguillon

Prediction of crack nucleation and propagation in porous ceramics using the phase-field approach

Keywords:

The relevance of the phase-field approach to simulate the fracture in porous ceramics has been investigated. For this purpose, the conditions for the crack initiation using the phase-field model have been compared to the theoretical predictions of the coupled criterion considering a pure V-notch singularity and a crack blunted by a cavity. For two types of ceramics (3YSZ and 8YSZ), it has been shown that the phase-field approach is able to simulate accurately the crack nucleation as predicted by the coupled stress-energy criterion. The nature of the regularization parameter l for the phase-field model has been discussed as function of the material and the local geometry where the crack initiates. Moreover, the apparent fracture toughness and the compressive fracture strength have been calculated for real porous YSZ ceramics. It has been found that the fracture properties of these complex 3D porous materials can be correctly predicted with the phase-field model. For specimen loaded under compression, it has been shown that the model is able to capture and explain the transition from a brittle behavior towards a diffuse damage when increasing the porosity.

Introduction

Porous ceramics are gaining a noteworthy interest as functional materials for many applications. They benefit from the ceramic properties, which are known to be highly stable in corrosive and oxidizing environments. Moreover, they offer high thermochemical and mechanical stability, low thermal conductivity and small dielectric constant [START_REF] Al-Naib | Introductory Chapter: A Brief Introduction to Porous Ceramic[END_REF][START_REF] Goto | A New Thick Film Coating Technology-Laser Chemical Vapor Deposition[END_REF][START_REF] Hammel | Processing and properties of advanced porous ceramics: An application based review[END_REF][START_REF] Ito | Heat-Resistant Coating Technology for Gas Turbines[END_REF][START_REF] Kawahara | Application of High-Temperature Corrosion-Resistant Ceramics and Coatings under Aggressive Corrosion Environment in Waste-To-Energy Boilers[END_REF][START_REF] Liu | General Introduction to Porous Materials[END_REF][START_REF] Liu | Applications of Porous Ceramics[END_REF][START_REF] Ohji | Testing and Evaluation of Mechanical Properties[END_REF]. Thanks to all these features, porous ceramics are used as key components in several technological devices with far-reaching economic and ecological implications. They are currently employed in biology, medical industry, electronic engineering, petroleum industry, metallurgy, environment protection and chemical engineering [START_REF] Goto | A New Thick Film Coating Technology-Laser Chemical Vapor Deposition[END_REF][START_REF] Hammel | Processing and properties of advanced porous ceramics: An application based review[END_REF][START_REF] Ito | Heat-Resistant Coating Technology for Gas Turbines[END_REF][START_REF] Kawahara | Application of High-Temperature Corrosion-Resistant Ceramics and Coatings under Aggressive Corrosion Environment in Waste-To-Energy Boilers[END_REF][START_REF] Liu | General Introduction to Porous Materials[END_REF]. In particular, porous ceramics have achieved an overwhelming success in the domain of energy. For example, they are used as barrier coatings to enhance the lifetime of many devices operated under severe conditions such as gas turbines. Porous ceramics have been also selected as material for the key components of several energy storage and power generation systems such as the solid oxide fuel cells and batteries [START_REF] Ghatee | Investigation of electrical and mechanical properties of 3YSZ/8YSZ composite electrolytes[END_REF][START_REF] Leonide | Long Term Operation of Rechargeable High Temperature Solid Oxide Batteries[END_REF][START_REF] Zhu | Advantages of intermediate temperature solid oxide fuel cells for tractionary applications[END_REF]. Yttria Stabilized Zirconia (YSZ) is a specific example, which exhibits good chemo-thermo-mechanical properties and high ionic conductivity. Nevertheless, despite all these advantages, the fracture of porous ceramic remains a problematic issue for all these technological applications.

It is worth noting that the mechanical behavior and especially the fracture of porous ceramics cannot be investigated based only on experimental studies because the analysis would require an unlimited number of mechanical characterizations. Indeed, the mechanical response of porous ceramics is complex and strongly dependent on the microstructure. For example, the apparent toughness of porous ceramic generally decreases with increasing porosity [START_REF] Samborski | Dynamic Fracture Toughness of Porous Ceramics: Rapid Communications of the American Ceramic Society[END_REF]. The crack onset and propagation is also a function of the local geometry of the pores and their distribution. Under specific conditions, the crack tip blunted by distant macropores can even lead to enhance the apparent toughness compared to the dense material [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF][START_REF] Deng | Reinforcement by crack-tip blunting in porous ceramics[END_REF]. For porous ceramics tested under compression, it has been shown that the fracture mechanism changes from a brittle behavior at low porosity towards a diffuse damage at high porosity [START_REF] Meille | Mechanical properties of porous ceramics in compression: On the transition between elastic, brittle, and cellular behavior[END_REF][START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. However, this transition has not been precisely explained yet. From this point of view, a numerical approach would be well adapted to better understand the precise role of the microstructure and to provide accurate prediction of the fracture of porous ceramics. Nevertheless, the complexity and the heterogeneity of the microstructure make the numerical treatment very challenging. For this reason, the simulation of the nucleation and propagation of cracks in such microstructure remains a salient subject of investigation. Therefore, there is still not a 3D model, which has been fully validated to predict precisely the fracture in porous ceramics.

To date, several works have been focused on solving the fracture problem in porous ceramics using different numerical approaches. On the one hand, it has been studied without simulating explicit cracks. In this case, phenomenological laws taking into account the loss of rigidity induced by the cracks nucleation and extension are used in order to reproduce the specimen macro-response during the loading [START_REF] Le Corre | Numerical modeling of the effective ductile damage of macroporous alumina[END_REF][START_REF] Smolin | Numerical simulation of mesomechanical behavior of porous brittle materials[END_REF]. Even if these methods enable the computation of the overall material behavior with good approximation, they are not based on physical models for fracture and their applications is generally limited to the studied case. On the other hand, methods based on the explicit simulation of cracks has been adopted in order to thoroughly study fracture initiation, propagation and branching. For this purpose, two main computational families have been developed in the frame of the Finite Element Method (FEM) depending on the sharp or diffuse description of the simulated crack. The first numerical approach is referred as the discrete crack model. In this case, the sharp crack is introduced as a discontinuity in the mesh and its propagation is simulated using various numerical methods. For instance, the extended finite element method (XFEM) [START_REF] Moes | A finite element method for crack growth without remeshing[END_REF], the phantom-node method, the element-erosion method, the generalized finite element method (GFEM) or re-meshing techniques [START_REF] Molnár | 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture[END_REF][START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Rabczuk | Computational Methods for Fracture in Brittle and Quasi-Brittle Solids: State-of-the-Art Review and Future Perspectives[END_REF][START_REF] Seleš | Crack propagation prediction in heterogeneous microstructure using an efficient phase-field algorithm[END_REF] have been used to simulate the fracture. However, the implementation of discrete crack models requires criteria for crack initiation and propagation. These criteria are taken from the Linear Elastic Fracture Mechanic (LEFM), which is well appropriate to describe the brittle behavior of ceramic materials. The extension of preexisting crack is governed by an energy balance based on the Griffith theory [START_REF] Gurtin | On the energy release rate in quasi-static elastic crack propagation[END_REF][START_REF] Gurtin | Thermodynamics and the griffith criterion for brittle fracture[END_REF]. For the fracture initiation, the coupled criterion proposed by D. Leguillon [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF] assumes that both the stress and energy criteria have to be fulfilled. It provides a rigorous framework for computing the crack nucleation [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Li | Finite element implementation of the coupled criterion for numerical simulations of crack initiation and propagation in brittle materials[END_REF]. However, the numerical implementation of the discrete crack models can be difficult since the fracture criterion must be evaluated on the whole structure in order to identify the sites for the crack nucleation. Furthermore, special algorithms are required to introduce discontinuities by modifying locally the shape functions or by re-meshing the path for the propagation. For all these reasons, discrete crack models are not well adapted to simulate complex cases like the branching and multiple cracks problem. Therefore, they cannot be applied to porous ceramics, which have complex and heterogeneous microstructures.

In order to overcome the limitations inherent to the discrete approach, several models have been proposed in the frame of the continuous fracture mechanics theory [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : I. Les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : II. Les modèles à gradient[END_REF][START_REF] Nguyen | Smoothing gradient damage model with evolving anisotropic nonlocal interactions tailored to low-order finite elements[END_REF][START_REF] Bui | A localized mass-field damage model with energy decomposition: Formulation and FE implementation[END_REF][START_REF] Vuong | Enhancement of the smoothing gradient damage model with alternative equivalent strain estimation for localization failure[END_REF][START_REF] Nguyen | Role of interfacial transition zone in phase field modeling of fracture in layered heterogeneous structures[END_REF]. All are based on a 'diffuse' or 'smeared' description of the discontinuity where the crack is implicitly modeled through a smooth scalar damage variable. In this context, Marigo and coauthors have proposed for the development of these models a rigorous theoretical framework using the variational approach and the gamma-convergence theory for the regularization of the free-discontinuity medium [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF][START_REF] Bourdin | The Variational Approach to Fracture[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : I. Les concepts fondamentaux[END_REF][START_REF] Pham | Approche variationnelle de l'endommagement : II. Les modèles à gradient[END_REF]. As a result, the fracture problem has been reformulated as the minimization of an energy functional without any additional criterion. On the basis of these developments, Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] have proposed a specific smeared crack model referred as the Phase-Field Method (PFM) for fracture mechanics. Thanks to the thermodynamic consistency of this model, a staggered scheme has been proposed for the numerical resolution [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF], allowing a simple implementation in FEM codes [START_REF] Seleš | Crack propagation prediction in heterogeneous microstructure using an efficient phase-field algorithm[END_REF][START_REF] Zhou | Phase field modeling of quasi-static and dynamic crack propagation: COMSOL implementation and case studies[END_REF][START_REF] Zhou | Phase field modelling of crack propagation, branching and coalescence in rocks[END_REF]. The crack extension is governed by an energy driving force that can be interpreted as a Griffith criterion. Several studies have shown that the model is able to simulate accurately the propagation, the branching and the multi-cracking π potential energy (J) Ψ density of elastic strain energy (J m -3 ) Ψ + density of strain energy related to tension contribution (J m -3 ) H 'history' loading field (J m -3 ) f thermodynamic driving force (J m -3 ) l length scale parameter for the phase-field approach (mm or µm) l 8YSZ length scale parameter for the 8YSZ material (mm or µm) l 3YSZ length scale parameter for the 3YSZ material (mm or µm) l Irwin

Irwin length (µm) l 8YSZ

Irwin

Irwin length for the 8YSZ material (mm or µm) l 3YSZ

Irwin

Irwin length for the 3YSZ material (mm or µm) l 0 initial crack length for the coupled criterion (mm or µm) l c 0 critical initial crack length for the coupled criterion (mm or µm) l c,8YSZ 0 critical initial crack length for the 8YSZ material (mm or µm) l c,3YSZ 0 critical initial crack length for the 3YSZ material (mm or µm) problem in complex geometry [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Seleš | Crack propagation prediction in heterogeneous microstructure using an efficient phase-field algorithm[END_REF]. For example, Nguyen et al. [START_REF] Nguyen | Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro CT imaging[END_REF] have successfully simulated the cracks evolution in a lightweight concrete microstructure with encapsulated pores. For all these reasons, the PFM seems suitable to study the fracture in porous ceramics. However, controversy remains on the PFM capacity to predict accurately the crack initiation.

On the one hand, Amor et al. [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF] have suggested that the PFM approach can predict the crack initiation by choosing an ad-hoc value for the length scale parameter l , which was initially introduced in the method to control the width of the smeared crack. Indeed, it has been proposed that this regularization length for the phase-field method can be related to the Irwin characteristic length linked to the material properties, i.e., strength and toughness. For this purpose, an expression of l has been derived considering a one-dimensional traction test [START_REF] Amor | Regularized formulation of the variational brittle fracture with unilateral contact: Numerical experiments[END_REF][START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF]. It can be noticed that some authors have successfully tested the relevance of this approach to simulate the crack nucleation for simplified geometries [START_REF] Nguyen | On the choice of parameters in the phase field method for simulating crack initiation with experimental validation[END_REF][START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]. Besides, Molnár et al. [START_REF] Molnár | Toughness or strength? Regularization in phase-field fracture explained by the coupled criterion[END_REF] have recently proposed a generalization of this expression for the 2D case still considering a homogeneous phase-field. They have found that the length scale parameter must be also a function of the ratio of the principal stresses suggesting a dependence to the local geometry where the crack nucleates. However, the generalization and the relevance of these expressions for complex 3D geometries without neglecting the gradient term in the phase-field is still questionable.

On the other hand, Kumar et al. [START_REF] Kumar | Revisiting nucleation in the phase-field approach to brittle fracture[END_REF] have recently criticized the above interpretation of l as a material constant. They have claimed that 'the fracture nucleation cannot be properly modelled by the phase-field formulation' since 'it is purely energetic' and does not take into account the material fracture strength. In this view, the regularization length is just a numerical 'parameter that is void of any further physical meaning'. To overcome this difficulty, they proposed to add an external driving force to the classical formulation depending on the material strength and supplementary numerical correction factors [START_REF] Kumar | Revisiting nucleation in the phase-field approach to brittle fracture[END_REF]. This literature review indicates that the interpretation of the regularization length and the use of the PFM to simulate the crack nucleation is still a subject of investigation. In this objective, the comparison of the PFM with the results given by a criterion based on the LEFM could be helpful to unravel the real nature of the length scale parameter l .

The aim of the article is to evaluate the relevance of the PFM approach to predict correctly the fracture in the complex 3D microstructures of porous ceramics. In this frame, the nature of the length scale parameter l is discussed depending on the material properties and geometry. For these objectives, the PFM simulations have been compared to the results of the coupled criterion obtained for stress singularity and stress concentration configurations (i.e. for a V-notched sample and a macro-crack blunted by a pore) [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF][START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavitiessize effects[END_REF]. The comparison has been conducted considering two types of zirconia stabilized with 3% and 8% molar of Yttria (3YSZ and 8YSZ) exhibiting different fracture properties. Once the capacity of the PFM model to detect accurately the crack initiation has been verified for these simplified geometries, the model has been used to compute the apparent fracture toughness and the compressive fracture strength of porous ceramics. The simulations have been compared to experimental data reported in [START_REF] Khajavi | Double Torsion testing of thin porous zirconia supports for energy applications: Toughness and slow crack growth assessment[END_REF] and [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF] as a function of the porosity for the model validation. Finally, the fracture mechanism for porous ceramics subjected to a compressive loading has been more precisely discussed.

Methodology

The objective of the present study is to evaluate the capability of the PFM approach to accurately predict the crack initiation in porous ceramic materials. For this purpose, PFM simulations have been conducted on different geometrical configurations and compared to the apparent fracture toughness given by the coupled criterion. First, Vnotched samples with different opening angles have been studied to address the case of crack nucleation from an ideal stress singularity [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF]. In order to mimic the fracture in porous ceramics, the crack nucleation from a stress concentration represented by a macro-crack blunted by a cavity [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF][START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavitiessize effects[END_REF] has been simulated. Moreover, the role of the material properties has been investigated considering two types of Zirconia Stabilized with 3% and 8% molar of Yttria (3YSZ and 8YSZ). Indeed, as reported in Table 1, the yttrium content in YSZ has a strong impact on the material fracture properties. For example, the strength and toughness of 3YSZ are roughly two and three times higher than 8YSZ, respectively [START_REF] Ghatee | Investigation of electrical and mechanical properties of 3YSZ/8YSZ composite electrolytes[END_REF]. The relevance of the PFM approach for predicting the propagation of a preexisting macro-crack in real porous ceramics has been assessed using 2D numerical microstructures. The results have been compared to the data reported in [START_REF] Khajavi | Double Torsion testing of thin porous zirconia supports for energy applications: Toughness and slow crack growth assessment[END_REF] related to an apparent fracture toughness measured on porous 3YSZ specimens. All the simulations have been performed with the plane strain elasticity assumption for a classical three-point bending symmetric test. In this condition, the crack is nucleated under a pure symmetric mode. Finally, the PFM has been used to simulate the crack initiation in real 3D porous ceramic microstructures for 8YSZ samples submitted to a compression loading. The results have been compared to the compressive fracture strength measured on the same microstructures reported in [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF].

Crack nucleation from a stress singularity: V-notched sample

The considered geometry of the V-notched sample with different opening angles 2β is illustrated in Fig. 1. The length and height of the beam are L = 8 mm and h = 2 mm, respectively, while the notch depth is a n = 0.4 mm. It is worth noting that the coupled criterion allows computing the apparent toughness K notch Ic for the notch while the critical load F c triggering the crack initiation is calculated with the PFM model. Therefore, the two approaches have been compared using the following expression [START_REF] Dunn | Fracture initiation at sharp notches: Correlation using critical stress intensities[END_REF]:

K notch Ic = ( 3F c L 2bh 2 ) h (1-λ) f ( a n h ) (1)
where b denotes the beam thickness taken to the unity for the 2D simulations. The exponent λ is the singularity order for a symmetric loading, which is obtained by solving the following equation [START_REF] Dunn | Fracture initiation at sharp notches: Correlation using critical stress intensities[END_REF]:

λsin(π -β) + sin(2λ(π -β) ) = 0 (2)
The dimensionless factor f ( an h ) in Eq. ( 1) is the polynomial function given in Eq. (3):

f ( a n h ) = c 1 a n h + c 2 a n h 2 + c 3 a n h 3 + c 4 a n h 4 + c 5 a n h 5 (3)
The coefficients c i are tabulated in [START_REF] Dunn | Fracture initiation at sharp notches: Correlation using critical stress intensities[END_REF] for a ratio an h ranging between 0.05 and 0.7. It can be noticed that the dimensions of the simulated specimen have been chosen to fulfill this condition on the ratio an h . Three notch angles 2 β ∈ {60 • , 90 large range of singularity orders with a crack initiation mainly controlled by the energy criterion for 2 β = 60 • and a non-negligible contribution of the stress criterion for 2 β = 120 • .

Crack nucleation from a stress concentration: crack blunted by a cavity

The simulated beam with a preexisting crack blunted by a cavity is shown in Fig. 2. Its geometry is the same as the previous V-notched case with a same crack length of a c = 0.4 mm. In this case, the fracture toughness K Ic and the applied critical loading F c are linked as follows:

K Ic = ( F c b ̅̅ ̅ h √ ) f ( a c h ) (4)
The dimensions of the studied geometry fulfil the condition ac h ∈ [0.05, 0.7] for which the dimensionless factor f ( ac h

) is expressed as [START_REF] Anderson | Fracture Mechanics: Fundamentals and Applications[END_REF]:

f ( a c h ) = 3 L h ̅̅̅ ac h √ 2 ( 1 + 2 ac h )( 1 -ac h )3 2 [ 1, 99 - a c h ( 1 - a c h ) { 2, 15 -3, 93 ( a c h ) + 2, 7 ( a c h ) 2 } ] (5)
Four porosities have been considered with a radius r ranging from 0.01 mm to 0.03 mm. These cavity sizes remain very small with respect to the crack length so that Eqs. ( 4) and ( 5) remain valid for assessing the apparent fracture toughness K app Ic . Moreover, the width δ c of the preexisting crack introduced in the mesh was small as compared to the pore size (i.e. δ c = 3 µm≪2r).

Crack nucleation and propagation in porous ceramics microstructures 2.3.1. Propagation of a preexisting macro-crack in 2D

The sample geometry to simulate the crack extension with the PFM approach in porous ceramic is displayed in Fig. 3a. The selected dimensions are L = 40 µm, h = 10 µm and a c = 3 µm. The ratio ac h is still comprised between 0.05 and 0.7 in such a way that Eqs. ( 4) and ( 5) can be used to express the toughness as a function of the applied loading. As shown in Fig. 3a-3c, three microstructures of different porosities ε ∈{10 %, 20 %, 30 %} have been simulated. These numerical microstructures have been generated using the Gaussian random field method [START_REF] Lanzini | Microstructural characterization of solid oxide fuel cell electrodes by image analysis technique[END_REF][START_REF] Moussaoui | Stochastic geometrical modeling of solid oxide cells electrodes validated on 3D reconstructions[END_REF].

Indeed, it has been shown that this method is able to emulate accurately the morphology of partially sintered ceramic microstructures, such as solid oxide fuel cell electrode materials [START_REF] Moussaoui | Stochastic geometrical modeling of solid oxide cells electrodes validated on 3D reconstructions[END_REF][START_REF] Moussaoui | 3D geometrical characterization and modelling of solid oxide cells electrodes microstructure by image analysis[END_REF]. The three microstructures have been generated by maintaining a constant correlation length for the solid YSZ phase [START_REF] Moussaoui | 3D geometrical characterization and modelling of solid oxide cells electrodes microstructure by image analysis[END_REF]. As result, the mean diameter for YSZ calculated by continuous Particle Size Distribution (PSD) [START_REF] Abrahamsen | A review of Gaussian random fields and correlation functions[END_REF], is approximatively a constant equal to 0.45 µm for all the microstructures simulated in the present study. On the contrary, the mean pore diameter decreases from 0.22 µm to 0.05 µm when densifying the ceramic from 30 % to 10 %. In all cases, the size of the porosity remains small compared to the simulated preexisting crack length. Besides, the crack width set to δ c = 0.1 µm is lower than the mean pore diameter for a porosity of 20 % and 30 %, whereas, for a porosity of 10 %, the two dimensions δ c and 2r are equivalent.

Crack nucleation in 3D

To study the crack nucleation in a real microstructure, the phase- field approach has been used to simulate our experiments reported in [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. In this previous work, the compressive fracture strength was measured as a function of the porosity ranging from 33% to 63%. The tested specimens were micro-pillars with a height of 130 µm made of porous 8YSZ. Due to a slight taper angle of 6 • , the pillars presented a section that evolves from 60 µm at the top to around 88 µm at the bottom. To measure accurately the porosity, 3D reconstructions were also obtained for all the investigated specimens [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF].

To predict the dependence of the compressive fracture strength with porosity, a multi-scale approach has been adopted coupling computation on the whole micro-pillar with simulation on the 3D porous ceramic microstructure (Fig. 4).

At the microstructure scale, the crack nucleation and propagation has been computed on sub-volumes, extracted from the whole reconstructions, using the PFM model. For intermediate porosities below 50%, a volume of 4 × 4 × 4 µm 3 has been simulated whereas a larger volume of 15 × 15 × 15 µm 3 has been considered for the highest porosity at 63%. These dimensions of the simulated domain were chosen to be representative of the heterogeneous microstructures. Indeed, the pore mean phase diameter is lower than 0.13 µm for the specimen at intermediate porosities (<50%) and equal to 0.99 µm at 63% [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. Therefore, even in this later case, several pores remain included along the edge of the computed volume (the characteristic size for the simulated volume is around fifteen times larger than the mean pore diameter). The solid phase of the reconstructed sub-volumes has been meshed using the software Avizo® considering tetrahedral elements. For the microstructures at intermediates porosities, the mesh were built using around 3.5 × 10 6 tetrahedrons (corresponding to around 17 × 10 6 degrees of freedom) while the volume at 63% contained 5.6 × 10 6 tetrahedrons (corresponding to around 26 × 10 6 degrees of freedom). A special attention was paid to the quality of the mesh by keeping elements with an aspect ratio below 10 and a dihedral angle above 10 • . For instance, a zoom of the generated mesh at a porosity of 63% is displayed in Fig. 4b. It can be seen that the surface of the solid phase inside the microstructure is smoothly meshed by the triangular elements avoiding numerical artifacts during the computation. The uniaxial compression has been simulated on the meshed reconstructions in displacement-controlled loading considering the properties of the dense 8YSZ (Table 1).

On the other hand, finite element simulations of the compression test at the micro-pillar scale were conducted assuming a homogeneous medium with a pure elastic behavior. For this modeling, a 2D approach including the full geometry of the tapered micro-pillar attached to a part of the pellet was considered using axisymmetric conditions (Fig. 4a). The size of the simulated pellet substrate was chosen in such a way that the mechanical response of the micro-pillar becomes independent of boundary conditions applied at the bottom of the substrate (i.e. the displacements are blocked in the axial direction). The simulations were performed using a Poisson's ratio of 0.31 and an effective Young's modulus taken from our experiments as a function of porosity [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. Because of the taper angle, it is worth noting that a stress gradient appears in the pillar during compression. For this reason, the cracks were observed during the experiments in the upper part of the pillar where the stress is the highest [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. The Volume Of Interest (VOI) related to the simulated domain at the microstructure scale was therefore taken at a short distance from the top of the pillar (i.e. 11 µm) where the damage is expected to occur (Fig. 4a). The difference in the displacements at the top and the bottom of the VOI surface was then extracted from the simulation at the macroscopic scale and applied on the 3D reconstruction simulated with the phase-field. In this approach, the compressive fracture strength corresponds to the critical load applied at the top of the micro-pillar identified when the damage is initiated in the simulated microstructure.

Models description

In this work, the coupled criterion has been used in order to validate the results obtained by the PFM model. The coupled criterion is thus briefly reminded hereafter before detailing the implemented PFM method.

Coupled criterion

The coupled criterion stipulates that the crack initiation is triggered when both the stress and energy criteria are fulfilled. This condition can be reached when the applied stress σ 0 reaches a critical value σ c 0 leading to the nucleation of a crack with an initial critical length l c 0 [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF]. When the crack initiation occurs with an initial length l 0 , the stress criterion requires that the local stress must exceed the material strength σ c along the expected crack path, as expressed in Eq. ( 6). Regarding the energy criterion, the incremental energy release rate G inc , which is the change in the potential energy W p due to the crack nucleation, must be higher than the material fracture toughness G c as shown in Eq. ( 7):

σ(x) ≥ σ c for 0 < x < l 0 (6) G inc (l 0 ) = - δW p l 0 ≥ G c (7)
The local stress σ(x) and the energy release rate G inc (l 0 ) introduced in the two last equations can be expressed as a function of the applied stress σ 0 using matched asymptotic expansions between the singular and the far fields [START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF]:

σ(x) = k(x)σ 0 for 0 < x < l 0 (8) G inc (l 0 ) = A(l 0 )σ 2 0 ( 9
)
where the dimensionless coefficients k and A depend on the crack length l 0 , the material properties and the local geometry where the crack initiates.

It is worth noting that the coefficient k decreases with increasing l 0 whereas A is an increasing function of l 0 . Therefore, combining Eqs. ( 6) and ( 8), the stress condition provides an upper bound for the nucleation since the applied loading σ 0 increases with l 0 (Fig. 5). The lower bound for nucleation is provided by the energy condition combining Eqs. [START_REF] Liu | Applications of Porous Ceramics[END_REF] and [START_REF] Ghatee | Investigation of electrical and mechanical properties of 3YSZ/8YSZ composite electrolytes[END_REF]. The coupled criterion corresponds to the configuration where the two conditions are fulfilled by minimizing the applied stress σ 0 denoted σ c 0 . In this condition, a unique solution for a critical crack length l c 0 is obtained, as shown in Fig. 5.

In the case of the V-notched sample submitted to a symmetric loading, the coupled criterion is formulated as follows [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF]:

K I ≥ K notch Ic with K notch Ic = ( G c Λ(β) ) 1-λ σ c 2λ-1 (10)
where K I is the generalized stress intensity factor and the term Λ(β) is a scaling coefficient depending on the notch opening angle. In this case, the critical length can be given through a precise analytical expression:

l c 0 = 1 Λ(β) G c σ c 2 = 1 Λ(β) Ê K Ic 2 σ c 2 ( 11
)
where Ê is the Young modulus for the plane strain condition ( Ê = E 1-ν 2 ). From Eq. ( 11), it can be noted that the critical length is proportional to the Irwin length (l Irwin = KIc 2 σc 2 ) through a factor depending on the notch opening and the Young modulus.

For a crack blunted by a cavity, it is worth mentioning that the matched asymptotic expansion for expressing σ(x) and G inc (l 0 ) has been conducted by taking into account the local geometry with the pore size [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF][START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavitiessize effects[END_REF]. Therefore, the apparent toughness for a blunted crack K app Ic is not only dependent on l 0 but also on the pore radius r [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF]:

K app Ic = ( G c D(r, l 0 ) ) 1/2 ( 12
)
where D is a function of r and l 0 . In this case, the critical length for the crack nucleation must depend on both the pore radius and the Irwin

length (i.e. l c 0 = f ( r, KIc 2 σc 2
)

). However, there is no analytical expression for this length, which must be calculated by solving a non-linear equation as detailed in [START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavitiessize effects[END_REF].

Phase-field method

The PFM approach described in Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF] has been adopted for this work. The model has been built considering quasi-static and isothermal conditions with the assumption of small strain. The variables of the problem are time dependent because the damage evolves during the loading. Nevertheless, the time step t is not 'real' and is introduced only through the loading increment in the simulation. Hereafter, the studied domain and the position are denoted V and x, respectively.

As mentioned in the introduction, the crack is modelled through a smooth scalar damage variable d(x) referred as the phase-field. Besides, a regularization length scale parameter l has been introduced to control the region of transition from the pristine state (d(x) = 0) to the fully broken state (d(x) = 1) (Fig. 6). In order to propose a specific expression for the smeared crack surface Γ l (d), depending on the damage variable, Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] have assumed that the evolution of the damage in onedimension can be approximated by an exponential function d(x) = e -|x| l . This hypothesis had allowed expressing Γ l (d) in 1D which has been directly extended to the three-dimension case as follows:

Γ l (d) = 1 2l ∫ V d 2 + l 2 ∇d⋅ ∇ddV with V⊂R δ , δ ∈ [1, 2, 3] (13)
It is worth noting that the minimization of the functional Γ l provides an approximation of the phase-field d(x), which is a solution of the associated Euler-Lagrange type equation:

d -l 2 Δd = 0 in V with ∇d .n = 0 on ∂V (14)
where n is the outward normal on ∂V considering the natural boundary conditions.

The PFM approach is mainly based on an energy balance. Indeed, the evolution of the potential energy during the loading π is the sum of three contributions, namely the rate of the stored elastic energy Ės , the rate of the work of the external forces Ṗc and the energy dissipation rate Ẇc due to crack propagation:

π = Ės + Ẇc -Ṗc ( 15
)
where the symbol Ẋ denotes the time derivative of X.

• The last term of this equation is classically expressed through the displacement field u and the applied force t on the edge of the domain:

Ṗc = ∫ ∂V t⋅ uds (16)
• The energy dissipation rate associated to the damage evolution can be expressed as the rate of the work W c required to create a diffuse crack in the volume:

W c = G c Γ l (d) and Ẇc = G c dΓ l (d) dt = G c ∫ V ( 1 l d -l Δd ) ḋdV (17)
where G c is the critical energy release rate of the material.

• The rate of stored energy is expressed by the integral of the density of elastic strain energy Ψ:

Ės = ∫ V ΨdV (18)
In the implemented model, it is assumed that the crack cannot propagate under compression. To fulfil this requirement, a spectral decomposition of the strain tensor ε is carried out and the density of strain energy is split into a sum of two contributions Ψ + and Ψ -, related to tension and compression, respectively:

Ψ = (g(d) + α)Ψ + + Ψ - (19)
where g(d) is a degradation function introduced to take into account the material softening during damage. Among several options [START_REF] Kuhn | On degradation functions in phase field fracture models[END_REF], g(d) is chosen, in this work, as a quadratic function with g(d) = (1 -d) 2 . It can be noticed that the numerical parameter α, which is assigned a very low positive value, is added to ensure the stability of the simulation when d tends to 1.

Thanks to the introduction of the degradation function, the stress tensor σ is dependent on the damage field d(x) as follows:

σ(ε, d) = ∂Ψ(ε, d) ∂ε (20)
Besides, the thermodynamic driving force f for crack extension is expressed in the model as the derivative of the elastic strain energy stored in the material with respect to the damage variable:

f = - ∂Ψ(ε, d) ∂d = 2(1 -d)Ψ + (21)
In this condition, the driving force for the crack propagation is controlled by the part of the strain energy related to the tension. Finally, the rate of the stored energy can be written as follows:

Ės = ∫ V [ σ : { 1 2 ( ∇ u + ∇ T u) } -f ḋ ] dV (22)
The governing equation to solve the problem is given by the rate of the potential energy Eq. ( 18) combined with Eqs. ( 16), ( 17) and [START_REF] Rabczuk | Computational Methods for Fracture in Brittle and Quasi-Brittle Solids: State-of-the-Art Review and Future Perspectives[END_REF]. For each increment of loading (or time), the minimization of this energy functional must be carried out on each variable d and u. In a staggered approach for the resolution, the minimization at a frozen damage for the phase-field (i.e. d(x) = cte) leads to solve a 'classical' elastic problem: ∇⋅σ(ε, d) = 0 for the domain V with σ⋅n = t on ∂V [START_REF] Seleš | Crack propagation prediction in heterogeneous microstructure using an efficient phase-field algorithm[END_REF] Then, the resolution at a given loading (i.e. u = cte) yields the equation for the phase-field computation when the crack propagates with ḋ > 0:

f = 2(1 -d)Ψ + = G c ( 1 l d -l Δd
)

with ∇d⋅n = 0 on ∂V [START_REF] Gurtin | On the energy release rate in quasi-static elastic crack propagation[END_REF] It can be noticed that this equation describes a classical energy criterion. Indeed, the crack is assumed to propagate when the driving force equivalent to an energy release rate in J⋅m -3 is equal to the density of energy created in the domain by the apparition of the cracked surface. Moreover, the formation of cracks is an irreversible process, so that the phase-field d(x) cannot decrease during the loading. To take into account this supplement constrain, Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] have proposed to introduce in Eq. ( 24) an 'history' loading field H (x, t) stipulating that Ψ + can only increase over the time. The partial differential Eq. ( 24) is thus rewritten as follows:

2(1 -d)H = G c ( d l -l Δd ) with H (x, t) = max s∈[0,t] Ψ + (u(x, t)) ( 25
)
The computation for the coupled criterion has been carried out with the FE code Modulef, while the phase-field approach has been implemented in the software Comsol Multiphysics®.

For the phase-field computations, it is worth noting that a special attention was paid to check the convergence of the simulations. More specifically, the number of iteration for the staggered resolutions was chosen so that the results are independent on the loading step. For this purpose, when the damage parameter exceeds 0.4, the loading step was divided by 100. Moreover, the independence of the result with the mesh size was also verified for all the simulations. It can be noticed that a very fine mesh was considered in the damage zone: for instance, more than 15 elements were included in the width 2l for the 2D simulations. Besides, all the computations have been carried out in such a way that the softening effect on the load-displacement curve due to a large damage zone remains negligible. Finally, the simulated tensile, shear and 3 point bending tests reported in Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] were used as benchmark to validate the implementation. For all the cases, a very good agreement has been found between the results from the Comsol Multiphysics ® and the data reported in Miehe et al. [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] validating the implementation.

Results and discussions

In this section, the results of the PFM simulations, conducted on the stress singularity and stress concentration configurations, are presented and compared to those performed using the coupled criterion. The capacity of the model to predict accurately the crack nucleation is evaluated and the nature of the length scale parameter l is discussed. Then, the numerical results obtained on 2D and 3D porous ceramic microstructures are compared to the experimental data reported in [START_REF] Khajavi | Double Torsion testing of thin porous zirconia supports for energy applications: Toughness and slow crack growth assessment[END_REF] and [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF] in order to investigate the relevance of the PFM model to simulate the fracture in porous ceramics.

V-notched 8YSZ sample

Preliminary computations have been conducted with the PFM model to simulate the crack pattern in the V-notched 8YSZ sample considering two length scale parameters (i.e.l 1 = 0.01 mm andl 2 = 0.04 mm). It can be noted that the classical undamaged condition (∇d .n = 0) was retained for the notch. As shown in Fig. 7a and7b, the crack initiates from the stress singularity and propagates in pure symmetric mode for both cases. Therefore, as already stated, the PFM is able to identify correctly the zone where the crack nucleates as well as its direction of propagation whatever the value of l [START_REF] Miehe | A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits[END_REF] (Fig. 7a and7b). Moreover, it can be seen that the crack becomes thicker when l is larger since this regularization parameter controls the width of the region of transition from the pristine material to the fully broken state. Fig. 7c shows that the force-displacement curve for l = 0.04 mm is slightly below the one calculated for l = 0.01 mm before the fracture. This result is due to the material softening effect during the loading, which is more pronounced when l is larger as previously mentioned.

As discussed in the introduction, the critical loading triggering the fracture was found to be strongly dependent on the choice of l . To investigate its role on the crack onset, the critical loading has been determined for several opening angles of the V-notched 8YSZ specimen. The critical loading F c was retrieved from the simulated forcedisplacement curves when the phase-field parameter d reaches ≈1 at the notch tip. The evolutions of the apparent facture toughness deduced from these computations are compared to the coupled criterion pre-Fig. 7. Visualization of the crack pattern initiated from a notch at an opening angle 2β of 90 • for: a) l 1 = 0.01 mm, b) l 2 = 0.04 mm. c) Force-Displacement curves considering the 8YSZ material at a V-notch angle 2β of 90 • forl 1 = 0.01 mm and l 2 0.04 mm. dictions in Fig. 8a. When the regularization parameter is increased, the apparent fracture toughness computed with the phase-field is increased especially at the highest opening angles. At the lowest angle (2 β = 60 • ), it can be noted in Fig. 8a that the apparent fracture toughness is almost independent of l . This behavior is explained since the singularity order at 2 β = 60 • is very close to the exponent 0.5 for a perfect crack (i.e. λ = 0.5122 for 2 β = 60 • ). In this condition, the crack initiation is almost governed by an energy criterion. As pointed out by Tanné at al. [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF], the phase-field method, which is based on a pure energetic approach for the crack propagation, is thus able to simulate the fracture initiation from a sharp singularity with a low sensitivity on the length scale parameter. On the contrary, at the highest opening angle, the apparent fracture toughness is strongly dependent on the choice of l for the simulations (Fig. 8a). Indeed, the low singularity order at 2 β = 120 • (i.e.λ = 0.6157) means that the fracture behavior is partially controlled by the stress criterion [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF][START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavitiessize effects[END_REF]. In this case, the prediction provided by the phase-field method, which does not consider the stress criterion, becomes a function of the regularization parameter. In other words, the choice of l is crucial for an accurate prediction of the crack onset when the fracture is partially controlled by the stress criterion. In these conditions, the regularization parameter could be seen as a pure numerical parameter that must be adapted for each configuration, as suggested by Kumar et al. [START_REF] Kumar | Revisiting nucleation in the phase-field approach to brittle fracture[END_REF].

Nevertheless, it has been possible to identify a single value of l denoted l 8YSZ in such a way that the apparent fracture toughness computed by the phase-field matches the coupled criterion prediction over the full span of investigated angles (Fig. 8a). Indeed, the PFM simulations obtained withl 8YSZ = 0.01 mm provides a precise prediction of the crack nucleation with a mean relative error lower than 2.4%. For this ad-hoc value of l , the PFM model is thus able to accurately reproduce the apparent toughness evolution when the fracture mechanism is changed from an energy controlled criterion to a coupled criterion. In other words, the length scale parameter allowing retrieving the coupled criterion predictions is found to be independent of the notch opening angle. The possibility to identify a single value of l 8YSZ suggests that this parameter could have a physical meaning and would contain the information to take into account the fracture properties involved in both criteria (i.e. toughness and strength). From this point of view, l 8YSZ could be related to a characteristic length of the material such as the Irwin length l Irwin or the critical length l c 0 involved in the coupled criterion (cf. Eq. ( 11)). This statement is in good agreement with Tanné et al. [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF] and suggests that, as for the uni-axial traction [START_REF] Borden | A phase-field description of dynamic brittle fracture[END_REF], the length scale parameter must be proportional to the Irwin length.

The critical loading leading to the fracture is plotted in Fig. 8b as a function of the opening angles 2β for the investigated length scale parameters. As expected, the critical force is almost independent of l at the lowest angles, whereas it decreases with increasing l at the highest angles (Fig. 8b-c). This dependence of the critical force (Fig. 8b) associated to the apparent fracture toughness (Fig. 8a) with the regularization parameter l at high opening angles was also observed in Tanné et al. [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]. This behavior can be interpreted through the relationship between the regularization parameter and the material characteristic length (l Irwin or l c 0 ). Indeed, if l is proportional to l Irwin or l c 0 , a large value of this parameter is related to a very low material strength. Since the fracture is partially governed at high opening angle by the strength criterion [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF], the critical force or the apparent fracture toughness must decrease with increasing l .

It can be noted that a decrease of the critical loading F c is observed when increasing 2β for the highest length scale parameter (i.e. l = 0.4 mm in Fig. 8b). In this condition, the material fracture strength must be very low due to the dependence of σ c with l . Therefore, the stress criterion is more easily fulfilled while the contribution of the energy criterion is enhanced. The force triggering the fracture is thus decreased with increasing 2β since the energy criterion is more easily fulfilled when the opening angle is large (note that K notch Ic is still an increasing function of 2β at l = 0.4 mm due to the evolution of the term h (1-λ) in Eq. ( 1)). This discussion reinforces the claim that, for a V-notched singularity, the regularization parameter in the phase-field model must be related to a material characteristic length l Irwin or l c 0 .

V-notched 3YSZ sample

The role of the material properties on the regularization parameter has been investigated by repeating the same study considering the 3YSZ instead of 8YSZ. The apparent fracture toughness calculated with the PFM model is plotted in Fig. 9 as a function of the V-notched opening angles for different values of l . The results are consistent with the previous analysis since the same evolutions of K notch Ic with 2β and l than the ones discussed for 8YSZ are observed. These evolutions are compared to the predictions given by the coupled criterion in Fig. 9. As for the 8YSZ ceramic, it has been possible to identify an ad-hoc length scale parameter l 3YSZ that allows retrieving the evolution of the theoretical fracture toughness with the opening angle. Indeed, the PFM results obtained forl 3YSZ = 0.02 mm match the coupled criterion predictions with a mean relative error on the apparent fracture toughness estimated to 1.9%. It is here worth noting that l 3YSZ is two times higher than l 8YSZ . This statement clearly shows that the regularization parameter changes with the simulated material. Moreover, the ratio l 3YSZ l 8YSZ = 2 for the phase-field approach is almost equal to the theoretical ratio

l 3YSZ Irwin l 8YSZ Irwin = l c,3YSZ 0 l c,8YSZ 0 = 1.8.
Therefore, this analysis confirms that when considering the crack onset from a singularity, the regularization parameter in the phase-field method is proportional to a material characteristic length l Irwin or l c 0 . This statement has far-reaching practical implications for the study of the V-notch configuration with the PFM model. Indeed, only a single study considering one material could be sufficient to determine the corresponding l for any other types of brittle materials.

Crack blunted by a cavity in 8YSZ and 3YSZ

To study the crack nucleation from a stress concentration, a crack blunted by a cavity has been simulated for the 8YSZ and 3YSZ ceramics. The relative toughness R = K app Ic /K Ic obtained with the PFM method as well as the theoretical predictions given by the coupled criterion are plotted as a function of the pore radii in Fig. 10. It can be noticed that the phase-field simulations have been carried out for various length scale parameters. As shown in Fig. 10, the dependence of the PFM results on the regularization parameter is less pronounced when the pore radius is lowered. As for the notch configuration, this evolution can be explained since the geometry tends towards a sharp crack entirely controlled by a pure energy criterion. On the contrary, the dependence on the apparent toughness calculated with the PFM becomes significant at the highest investigated pore radius for which the contribution of the stress criterion to the fracture is non-negligible. In this case, the apparent fracture toughness is lowered with increasing l . This evolution can be interpreted through a very low facture strength associated to the large length scale parameter. This explanation would mean that the regularization parameter can be still related to a material characteristic length for a stress concentration [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]. As illustrated in Fig. 10a for 8YSZ, this dependence can lead to an incoherent behavior for l = 0.04 mm with a weakening effect with increasing the pore radius.

For both materials, it has been possible to identify two specific values for the regularization parameters l 8YSZ and l 3YSZ to fit accurately the predictions given by the coupled criterion (Fig. 10). Indeed, the mean relative error the relative fracture toughness simulated with l 8YSZ = 0.005 mm and l 3YSZ = 0.007 mm is equal to 1.1% and 1.2% for 8YSZ and 3YSZ respectively. Therefore, for a given value of l , the PFM model is able to capture the expected strengthening effect due to the crack blunted by a cavity [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF]. Nevertheless, it can be noted that the length scale parameters l 8YSZ and l 3YSZ identified for the V-notch singularity are not the same than the ones obtained for the blunted crack (for the same specimen submitted to a three-point bending test). For instance, the regularization parameter identified for the notch (l 8YSZ = 0.01 mm) is significantly higher than the one for the blunted crack (l 8YSZ = 0.005 mm). When the length scale parameter for the notch is taken to l 8YSZ = 0.005 mm, the phase-field results overestimate the coupled criterion predictions. In this case, the discrepancy between the critical forces provided by the coupled criterion and the phase-field reaches almost 20% for a notch opening angle of 2β = 120 • . Such a discrepancy cannot be explained by the numerical uncertainties for the phase-field simulations, which have been checked to be very limited in our case. This result would thus suggests that the regularization parameter could be also a function of the type of local geometry where the crack nucleates.

It can be noticed that the ratio l 3YSZ l 8YSZ = 1.4 obtained with the phasefield method for the blunted crack is different from the Irwin one l 3YSZ Irwin l 8YSZ Irwin = 1.8. However, it is almost equal to the ratio of the critical crack length in the coupled criterion l c,3YSZ 0 l c,8YSZ 0 = 1.5. This statement means that the length scale parameter l for the phase-field method is related to the critical length l c 0 for the coupled criterion. For a V-notch singularity, the Fig. 9. Dependence of the apparent fracture toughness on the opening angle for the 3YSZ material considering different values of the length scale parameter l . l 3YSZ refers to the length scale parameter that provides close agreement between the coupled criterion and the PFM results over the whole considered range of notch angles. 0 is equivalent to the ratio given by the Irwin characteristic length as mentioned in the previous section. However, for a crack blunted by a cavity, the critical crack length l c 0 is still a function of l Irwin but also takes into account the local geometrical characteristics of the stress concentration [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF][START_REF] Leguillon | Prediction of crack initiation at blunt notches and cavitiessize effects[END_REF] (cf. Section 3.1). This discussion reinforces the proposition that the length scale parameter must be dependent on the type of local geometry where the crack nucleates for a stress concentration. It can be mentioned that the possibility could be further discussed by comparing the coupled criterion prediction with alternative phase-field methods [START_REF] Nguyen | Smoothing gradient damage model with evolving anisotropic nonlocal interactions tailored to low-order finite elements[END_REF][START_REF] Bui | A localized mass-field damage model with energy decomposition: Formulation and FE implementation[END_REF][START_REF] Vuong | Enhancement of the smoothing gradient damage model with alternative equivalent strain estimation for localization failure[END_REF][START_REF] Nguyen | Role of interfacial transition zone in phase field modeling of fracture in layered heterogeneous structures[END_REF]. Nevertheless, even if additional studies are still needed, it can be concluded from this analysis that the PFM model can be applied to predict the crack initiation in porous ceramic with a length scale parameter that have to be identified for each type of microstructures.

Crack nucleation in porous microstructures

When the length scale parameter l is correctly chosen, the above case studies have shown the capacity of the PFM model to accurately predict the crack nucleation on ideal geometries (i.e. V-notch and the blunted crack). In this section, the relevance of the PFM model to simulate the fracture in porous media with complex microstructures is evaluated.

Apparent fracture toughness of porous ceramics

To study the evolution of the apparent fracture toughness with the porosity, 2D porous 3YSZ microstructures have been generated with the random field method considering a constant correlation length for the solid phase (cf. Section 2.3 and Fig. 3). As a result, the geometrical features of the simulated microstructures are preserved except modifications affecting mainly the porosity (i.e. volume fraction and pore mean diameter). Therefore, based on the previous statements, it is anticipated that a unique length scale parameter for the phase-field model is sufficient to predict the evolution of the apparent toughness with the porosity. Two methods have been used to determine K app Ic from the phase-field simulations. In the first one (denoted 'method 1 ′ thereafter), the toughness was deduced when the fracture initiates at the tip of the preexisting macro-crack (with δ c = 0.1 µm) blunted or not by a pore (cf. Fig. 11). In the second method (called 'method 2 ′ ), the toughness was assessed after propagation, when the fracture is reinitiated from the first pore reached by the crack. The reason for using two methods is to provide an estimation of the toughness scattering induced by heterogeneities in the simulated microstructures.

For each porosity, two independent realizations (i.e. microstructures) exhibiting the same geometrical statistical properties have been generated using the Gaussian random field method. As an illustration for a porosity of 20%, it can be seen in Fig. 12 that the PSD of the two microstructures shown in Fig. 11a-b are nearly perfectly superimposed. The simulations of these two microstructures at the same porosity level should allow estimating the impact of the local spatial variations on the toughness predictions. It is worth noting that, in spite of these statistical variations on the local geometry, the same type of morphology is preserved since the same method and input parameters have been used to generate the microstructure. Finally, as pointed out in Section 2.3.1, the synthetic microstructures generated with the Gaussian random field method provide a good approximation of real porous ceramic microstructures produced by powder sintering [START_REF] Moussaoui | Stochastic geometrical modeling of solid oxide cells electrodes validated on 3D reconstructions[END_REF]. Consequently, the experimental toughnesses reported in [START_REF] Khajavi | Double Torsion testing of thin porous zirconia supports for energy applications: Toughness and slow crack growth assessment[END_REF], which were measured using double-torsion testing on 3YSZ porous membranes produced by tapecasting, have been selected for comparison with the phase-field simulations.

All the simulations have been performed with a single length scale parameter l 3YSZ = 0.1 µm that remains small compared to the correlation length used to generate the microstructures (cf. Section 2.3.1). The apparent fracture toughnesses calculated on the first synthetic microstructure are compared to the experimental data in Fig. 13a. Independently of the method used to assess the toughness, the phasefield results are consistent with the measurements. The offset between the simulations and the experimental curve can be explained by the assumptions considered for the modeling, which are mainly (i) the differences between the synthetic microstructures and the real ones, and (ii) the hypothesis of 2D simulations for the phase-field model. Moreover, it is worth noting that the slight disagreement between the experimental data and the simulations can be further reduced by increasing l 3YSZ . For example, it has been found that a regularization parameter taken at l 3YSZ = 0.8 µm allows fitting almost perfectly the experimental data (indeed, at a porosity of 10 %, the computed fracture toughness with the method 1 falls to K app Ic = 3.8 MPa.

̅̅̅̅ m √

). However, a tradeoff has to be found for the length scale parameter. On the one hand, this parameter must be adapted to describe accurately the initiation step in order to predict the apparent fracture toughness. On the other hand, a sufficiently low value for l must be ensured to describe correctly the subsequent crack propagation in the fine microstructure. From a practical point of view, if necessary, two length scale parameters may be considered for the simulations. The first one could be adapted to compute precisely the fracture initiation, while the second one could be chosen small compared to the size of the microstructure in order to describe the crack propagation in the porous ceramic.

Nevertheless, with a single length scale parameter, the model captures correctly the decrease of K app Ic with increasing the porosity. This result confirms the relevance of the phase-field method to predict the fracture toughness of porous brittle materials after calibration of the regularization parameter for a given type of microstructure. This supports the claim that this parameter contains material fracture properties so that the PFM model can be used to quantify the conditions for the propagation of a pre-existing crack in porous ceramics.

The comparison of the fracture toughness predictions as a function of the porosity for the two microstructures (displayed in Fig. 12) is shown in Fig. 13b. For the sake of clarity, only the results with the first method have been reported. The results fall within the same range for the two simulated microstructures. The slight difference between the two curves is due to local variations in the geometry at the preexisting macro-crack tip. Therefore, by keeping the same length scale parameter for the phasefield model, it is possible to estimate correctly the fracture toughness for two different microstructures exhibiting similar statistical features for the morphology. In other words, the fracture properties can be determined with a single regularization parameter for one type of porous media (characterized by the same morphological features related to the manufacturing process or the method of synthetic microstructure generation). This statement is consistent with the previous discussion on the dependence of the regularization parameter on the type of the local geometry where the crack initiates. For a complex microstructure, l would thus depend on the characteristics of the solid and pore phases (i. e. shape and size of pores and of solid ligaments between the pores, their statistical distribution, etc. [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF]).

In Fig. 13a, it can be noticed that the first method to compute the apparent fracture toughness leads to slightly overestimate K app Ic compared to the second one at ε = 20 % and 30 % and underestimate K app Ic at 10 %. These results can be interpreted by inspection of the cartographies displaying the crack path within the microstructure in Fig. 11. For instance, at ε = 20 %, it can be noted that the preexisting macro-crack ends in a large pore (Fig. 11a). Due to the blunting effect induced by this large cavity, the apparent toughness computed with the first method is thus enhanced. On the contrary, after propagation, the crack reaches a smaller pore. Therefore, the apparent toughness for the crack re-initiation from this smaller cavity is lowered knowing that the other conditions affecting the fracture remains roughly similar. Indeed, the direction of propagation remains in a quasi-pure opening mode and the distance between the pores is almost constant (Fig. 11a) [START_REF] Leguillon | Fracture of porous materials -Influence of the pore size[END_REF]. At the lowest porosity fraction (ε = 10 %), the initial macro-crack tip is located in the 3YSZ solid phase (Fig. 11c). Therefore, in this particular case, the crack tip is no longer blunted by a cavity explaining the lower value of K app Ic obtained with the first method with respect to the second one. Moreover, it can be noticed that the direction of propagation remains similar (only the distance between the initial macro-crack tip and the first pore and between the first and second pores is changed).

Compressive fracture strength of porous ceramics

To study the crack nucleation in pristine porous ceramics, the compressive fracture strength measured in [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF] as a function of porosity for 8YSZ has been calculated with the PFM model. The simulations have been conducted on 3D porous ceramic microstructures according to the methodology detailed in Section 2.3.2. Two length scale parameters of l = 0.4 µm and l = 0.8 µm have been considered for the computations.

The calculated compressive fracture strength is compared to the experimental data as a function of porosity in Fig. 14. It can be seen that the numerical predictions are in very good agreement with the measurements. Especially at l = 0.8 µm, the model captures almost perfectly the decease of the compressive fracture strength with the porosity. This statement confirms that, with a single value of the regularization parameter l , the phase-field method is also able to predict accurately the crack nucleation in partially sintered ceramics. Therefore, for a given type of microstructure for which the regularization parameter has been identified, the phase-field approach can be seen as a predictive tool to study the fracture behavior of uncracked porous ceramic. In other words, as already discussed on ideal geometries, the length scale parameter must contain the material properties so that the phase-field model can be used to predict the fracture in complex porous ceramic microstructures.

It is worth noting that a transition in the fracture mechanism with the porosity fraction was detected in the micro-compression tests [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. Indeed, at low porosity, it was found that the fracture is controlled by a Fig. 13. Experimental and simulated dependences of the apparent fracture toughness on porosity using a) the two computation methods for realization 1, b) "method 1" for the two realizations. pure brittle behavior whereas a diffuse damage with the generation of stable micro-cracks in the microstructure was observed at high porosity. This behavior is well retrieved by the PFM model and can be analyzed thanks to the simulated results obtained on the 3D volumes. For instance, as shown in Fig. 15a-15c for a porosity of 33%, the apparition of the first crack in the microstructure is almost instantaneously followed by a complete damage spreading in the whole microstructure. Indeed, the first cracks are detected in the volume for an applied force on the micro-pillar of 6.04 N (Fig. 15a). Then, for a very small increment in charge (less than 25%), the damaged variable takes a value of d = 1 in the whole simulated domain indicating the total fracture of the specimen as shown in Fig. 15c. Therefore, as soon as the facture is initiated in the volume, it leads to the complete collapse of the material as it was observed during the experiments. It can be also noticed that the first micro-cracks are roughly aligned and parallel to the direction of the applied solicitation as detected during the tests [START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF]. On the contrary, for the highest porosity, the simulations have revealed the formation of a network of stable micro-cracks distributed in the microstructure during the loading (Fig. 15b-15d). As shown in Fig. 15b, the first micro-cracks appear for an applied force of 0.09 N. Then, the density of cracks is progressively increased with increasing the loading. The repartition of micro-cracks in is for example shown in Fig. 15d for an applied force of 0.16 N. The total collapse arises only when the load is increased up to 0.17 N. These results, which are in very good agreement with the experiments, constitutes a further proof of the relevance of the phasefield approach to calculate accurately the fracture initiation in complex porous ceramics.

The modeling results can be used to analyze this transition in the fracture mechanism as a function of the porosity. For this purpose, the density of elastic strain energy ξ stored in the microstructure has been calculated just before the first crack nucleates in the volume:

ξ = 1 V YSZ ∫ VYSZ ΨdV (26)
where V YSZ is the volume of the solid phase. This density of energy is plotted in Fig. 16 as a function of the porosity. It can be noticed that the strain energy strongly decreases with increasing the porosity. The high stored energy at low porosity should thus constitute the driving force to propagate the first nucleated micro-cracks to the whole microstructure explaining the brittle behavior of the ceramic. At the opposite, for the highly porous material, the stored energy becomes insufficient to ensure the propagation to the whole specimen. In this condition, the fracture is based on a diffuse damage as shown by the modeling and experimental results.

Moreover, the mean value of the local maximum principal stress, which was taken in the ligament where the first crack appears, has been also retrieved from the simulated data just before the fracture (i.e. the loading step corresponding to the time when the stress is the highest in the ligament before relaxation due to the fracture). It can be noticed that this principal stress is in traction. Indeed, due to the complexity of the microstructure loaded under compression, some part of the ceramics are submitted to a bending moment leading to local tensile stress. It is worth noting that this local tensile stress triggering the fracture is roughly equal to 4 GPa whatever the investigated porosities (i.e. 33 %-63 %). This high value could be consistent with the apparent fracture strength of the 8YSZ ceramic when considering the very small volume of the ligament submitted to the tensile loading [START_REF] Doitrand | Revisiting the strength of micronscale ceramic platelets[END_REF].

As for the ideal 2D geometries, the above discussion shows that the phase-field model is able to capture a change in the fracture criterion thanks to an accurate identification of the regularization parameter. For a porous ceramic loaded under compression, it enables us to explain the observed transition in the fracture mechanism from a brittle behavior towards a diffuse damage. From all these results, the phase-field model appears as a relevant tool to quantify the crack nucleation in porous ceramics but also to analyze the underlying mechanisms controlling the fracture.

Conclusion

The relevance of the phase-field model for simulating the crack nucleation and propagation in porous ceramics has been evaluated. For this purpose, the PFM results have been compared to the predictions of the coupled criterion based on linear fracture mechanics. The study was conducted on two types of ceramics, whose fracture properties are significantly different (i.e. 3YSZ and 8YSZ). A pure singularity induced by a V-notch and a stress concentration related to a crack blunted by a cavity have been studied considering a three point bending test.

In the V-notch case, considering a unique length scale parameter for the phase-field model l , it has been shown that the PFM method is able to retrieve the coupled criterion predictions as a function of the notchopening angle. Moreover, the ratio of the regularization parameter for the two ceramics is equal to the ratio for both the Irwin length and the critical length for the coupled criterion ( l 3YSZ l 8YSZ ≈ ). Therefore, the length scale parameter depends on the material fracture properties as expected. For an ad-hoc value of l , the PFM method is thus able to predict correctly the facture initiation at singularities even when the crack nucleation is controlled by a coupled energy and stress criterion. For a crack blunted by a cavity, the PFM model is also able to reproduce the theoretical dependence of the apparent fracture tough- ness on the pore size. In this case, the ratio of the selected length scale parameters for the two materials is no longer equal to the ratio for the Irwin length but remains identical to the one related to the coupled criterion ( l 3YSZ l 8YSZ ≈ l c,3YSZ

0 l c,8YSZ 0 ∕ = l 3YSZ Irwin l 8YSZ

Irwin

). Moreover, the length scale parameter for the V-notch was found significantly different from the one identified for the blunted crack. Therefore, this result suggest that the length scale parameter is also dependent on the type of local geometry or microstructure where the crack initiates. In this context, it should be determined using theoretical or experimental data. This analysis conducted with the coupled criterion and the PFM on ideal geometries indicates that the phase-field should be an appropriate method to predict the crack initiation in complex porous ceramic microstructures. To confirm this claim, the relevance of the phase-field method to simulate the crack nucleation and propagation in representative porous ceramic microstructures has been investigated. For this purpose, the conditions triggering the propagation of a pre-existing crack in porous 3YSZ have been computed with the phase-field method on 2D synthetic microstructures. As proposed, the apparent fracture toughness was found to be in good agreement with experimental data over the whole porosity range using a unique length scale parameter. Moreover, the conditions of crack re-initiation from porosity during the propagation has been discussed as a function of local microstructural parameters (i.e. pore size, length of the ligaments between the pores and direction of propagation). These results confirm that the phase-field method is able to predict accurately the fracture toughness of porous brittle materials after calibration of the regularization parameter.

The fracture initiation in uncracked porous ceramics have been studied on real 3D microstructures. In that objective, the compressive strength has been computed with the PFM model for 8YSZ with porosities ranging from 33% to 63%. For a single length scale parameter, it has been found that the compressive strengths computed as a function of the porosity are in very good agreement with the experimental data. Besides, it has been shown that the model predicts the transition from a brittle behavior at low porosity towards a diffuse damage at high porosity. This change in the fracture mechanism has been analyzed with the model. It has been shown that this transition is related to the elastic strain energy stored in the microstructure. At low porosity, the excess of strain energy allows the unstable crack propagations in the whole specimen. Conversely, at high porosity, the strain energy is insufficient to propagate the nucleated micro-cracks to the whole microstructure. All these results show that the PFM method is a relevant tool to predict accurately the fracture of porous ceramics.

Fig. 1 .

 1 Fig. 1. Three-point-bending test on a V-notched sample.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Three-point-bending test on a sample with a crack blunted by a cavity.

Fig. 4 .

 4 Fig. 4. Mesh of the: a) tapered pillar attached to the pellet with axisymmetric conditions. A zoom on the pillar is shown together with the location of the Volume Of Interest (VOI) for the simulation at the microstructure scale, b) 3D microstructure at a porosity of 63 %.

Fig.Fig. 6 .

 6 Fig. Applied loading as a function of the initial crack length for the coupled criterion obtained using the stress (in blue) and the energy (in purple) criteria. The coupled criterion is chosen among the admissible solutions such l c 0 minimizes the applied loading σ c 0 (in red). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 8 .

 8 Fig. 8. Dependence of the mechanical response of the V-notched 8YSZ sample on the opening angle considering different values of the length scale parameter l : a) apparent fracture toughness, b) critical force triggering crack initiation. l 8YSZ refers to the length scale parameter that provides close agreement between the coupled criterion and the PFM results over the whole considered range of notch angles, c) force-displacement curves at a V-notch angle 2β of 120 • .

Fig. 10 .

 10 Fig. 10. The dependence of the relative toughness on the cavity radius considering different values of the length scale parameter l for: a) 8YSZ material, b) 3YSZ material. l 8YSZ and l 3YSZ refer to the length scale parameters that provide close agreement between the coupled criterion and the PFM results over the whole considered range of cavity radii.

Fig. 11 .

 11 Fig. 11. Cartographies displaying the crack path within the porous microstructures: a) first realization at ε = 20 %, b) second realization at ε = 20 %, c) first realization at ε = 10 %. (l 3YSZ = 0.1 µm).

Fig. 12 .

 12 Fig. 12. Comparison between two generated synthetic microstructures with ε = 20 %: a) visualization of the two realizations of the microstructure, b) Solid phase size distributions computed on the two realizations.

Fig.

  Fig.Comparison between the calculated compressive fracture strength and the experimental[START_REF] Abaza | Fracture properties of porous yttria-stabilized zirconia under microcompression testing[END_REF] as a function of porosity.

Fig. 15 .

 15 Fig. 15. Visualization of the cracks created in porous microstructures submitted to compression: First micro-cracks created in the microstructure for a porosity of a) 33 % and b) 63 %. Evolution of the damage variable after a small increment of charge at c) 33 % and d) 63 % (l = 0.4 µm).

Fig. 16 .

 16 Fig. 16. Evolution of the density of elastic strain energy as a function of porosity (l = 0.8 µm).

  • , 120 • } have been investigated to sweep a

	Table 1
	Mechanical properties of 3YSZ and 8YSZ.
	Material	E [GPa]	ν	KIc[MPa.	√	̅̅̅̅ m ]	σc[MPa]	lIrwin =	( KIc σc	) 2	[µm]	Reference
	3YSZ	214	0.31	5.1			583	76.5				[59]
	8YSZ	216	0.31	1.61			245	43.2				[60,61]

Acknowledgements

This project has received funding from the Fuel Cells and Hydrogen 2 Joint Undertaking (JU) under grant agreement n • 825027 (AD ASTRA project). The JU receives support from the European Union's Horizon 2020 research and innovation program and Denmark, France, Italy, Spain, Poland, Netherlands, Greece, Finland, Estonia, Germany, United Kingdom, Switzerland. The work has also been partially supported by the Carnot Institutes (CASSIOPEE project) and Genvia company.

CRediT authorship contribution statement

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.