Kajetan Wojtacki

Pierre-Guy Vincent

Pierre Suquet
email: suquet@lma.cnrs-mrs.fr

Hervé Moulinec

Guylaine Boittin

A micromechanical model for the secondary creep of elasto-viscoplastic porous materials with two rate-sensitivity exponents: application to a mixed oxide fuel

Keywords: porous media, viscoplasticity, FFT method, homogenization, mathematical morphology, microstructures

est destinée au dépôt

A micromechanical model for the secondary creep of elasto-viscoplastic porous materials with two rate-sensitivity exponents: Application to a mixed oxide fuel

Introduction

Mixed oxide (MOX) fuel is a nuclear fuel that consists of plutonium oxide mixed with natural, depleted or reprocessed uranium oxide. MOX fuel has been used in French Pressurised Water Reactors since 1987 [START_REF] Oudinet | Characterization of plutonium distribution in MIMAS MOX by image analysis[END_REF]. This material exhibits a very specific microstructure which evolves significantly during irradiation [START_REF] Noirot | Detailed characterisations of high burn-up structures in oxide fuels[END_REF]. As an example, the accumulation of rare gases like Krypton or Xenon generated by the fission process, triggers the onset of irradiation bubbles within the material. This voiding in the bulk of the material has a direct influence on its physical properties and consequently on the loading applied by the fuel on the cladding. Therefore, in order to assess the safety of nuclear power plants, especially under accident conditions such as a Reactivity Initiated Accident (RIA), the impact of those microstructural differences on the overall thermo-mechanical behavior of the fuel has to be carefully studied [START_REF] Koo | A comparative analysis of UO2 and MOX fuel behavior under reactivity initiated accident[END_REF][START_REF] Schmitz | High burnup effects on fuel behaviour under accident conditions: the tests CABRI REP-Na[END_REF][START_REF] Sasajima | Behavior of irradiated ATR/MOX fuel under reactivity initiated accident conditions[END_REF][START_REF] Fuketa | Transient response of LWR fuels (RIA)[END_REF].

Our main interest here is in the creep properties of this fuel. A special attention is paid to creep tests in which the material is subjected to a constant stress, for fuel pellets mostly in compression and at high temperature. The response of the material is first an instantaneous elastic deformation, followed by an increase of the strain in a progressive and nonlinear way (primary creep) until the strain rate stabilizes (secondary or steady-state creep) [START_REF] Routbort | Compressive creep of mixed-oxide fuel pellets[END_REF]. Only the regime of secondary creep is addressed here.

The purpose of this work is to estimate both numerically and analytically the effect of clusters of voids on the strain rate of an irradiated MOX fuel at high temperature, 1573 K ≤ T ≤ 2473 K, and high stress level, 10 MPa ≤ Σ ≤ 100 MPa, in the range of conditions which could be encountered in a fuel during a postulated RIA [START_REF] Suzuki | Thermal stress analysis of high-burnup LWR fuel pellet pulseirradiated in Reactivity-Initiated Accident conditions[END_REF].

The overall creep behavior of the material is simulated by means of computational periodic homogenization. Three-dimensional full-field simulations are performed by a method based on Fast Fourier Transforms (FFT) using a home code (CraFT) freely available on the web1 . This mesh-free approach, first proposed by [START_REF] Moulinec | A fast numerical method for computing the linear and the nonlinear mechanical properties of composites[END_REF], is especially attractive for complex heterogeneous microstructures with non-linear phases. Full-field simulations require the reconstruction of artificial microstructures, since 3D images of irradiated nuclear fuels are not available. In comparison to the real microstructure, the artificial realizations have to match appropriate morphological descriptors [START_REF] Jeulin | Random texture models for material structures[END_REF], which have been recognized in a preliminary analysis as containing the most relevant information from the point of view of the considered physical phenomena. The geometry of these microstructures is simplified by considering only spherical pores located within randomly distributed spherical clusters. The matrix surrounding the pores is considered as isotropic and its viscoplastic behavior is described by a strain-rate potential including two creep mechanisms, diffusion creep and dislocation creep. The constitutive relations for the matrix are derived from the phenomenological strain rate identified for MOX fuel in a creep test by [START_REF] Slagle | Fuel transient deformation[END_REF]. This potential is described by two power-law functions with two different rate-sensitivity exponents.

The analytical model makes use of a variational approach and provides a lower bound on the effective strain-rate potential [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF][START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF][START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF][START_REF] Ponte Castañeda | Nonlinear composites[END_REF].

The paper is organized as follows. The microstructure of the material under study is analyzed in Section 2. A few simplifications relative to its morphology and the constitutive relations of the matrix material are introduced in Section 3 together with the procedure followed to generate the microstructures used in the full-field simulations. Section 4 gives a brief introduction to the FFT method and discusses the spatial resolution adopted in the full-field simulations. The analytical model is derived in Section 5. Finally, the accuracy of the model is assessed in Section 6 by comparing its prediction with full-field simulations.

Microstructure of MIMAS MOX fuel

As a result of their production process, MOX fuels are multi-phase materials composed of a matrix of uranium dioxide containing plutonium-rich agglomerates. During irradiation, the microstructure of this heterogeneous material evolves and these agglomerates become porous due to the accumulation of rare gases and to the apparition of irradiation bubbles.

In order to relate the overall viscoplastic behavior of this porous material to its microstructure, several morphological informations about the constituents of the material are extracted from data and micrographs from the literature. Since our interest is in the generation of numerical realizations mimicking the actual material, the morphological descriptors used in the present study to analyze the phases are (following [START_REF] Jeulin | Random texture models for material structures[END_REF]) their volume/surface fraction, their size (cluster or pore diameter, characteristic length) and their spatial distribution. The geometrical description of the MOX microstructure is performed below in two steps. First the agglomerates (or clusters) are analyzed and second the pores (or irradiation bubbles) are studied.

Pu-rich clusters

The fabrication process has a direct impact on the morphology of the fuel. The present study concerns MIMAS MOX fuels. The MIMAS process (MIcronized MASter blend) is widely used for MOX fuel fabrication [START_REF] Oudinet | Characterization of plutonium distribution in MIMAS MOX by image analysis[END_REF]. It consists in a micronisation of a mixture with 25-30% plutonium dioxide (master-blend mixture), followed by dilution with natural or depleted uranium dioxide to reach the desired plutonium content (more details on the fabrication routes are given in [START_REF] White | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF] and [START_REF] Fisher | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF]). The microstructure of the MOX fuel obtained by the MIMAS process depends also on the characteristics of the uranium dioxide powders, which can be produced by AUC (Ammonium Uranate Carbonyl) or ADU (Ammonium DiUranate) routes. The manufacturing process returns a material with plutonium-rich agglomerates. The MIMAS-AUC has a classical duplex structure with about 25 vol% Pu agglomerates in a UO 2 /PuO 2 lattice, whilst the MIMAS-ADU has three different phases, the UO 2 rich matrix, the Pu-rich agglomerates, and a coating phase with intermediate Pu content around the particles of UO 2 powder [START_REF] White | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF].

The maximal size of Pu-rich spots reported in [START_REF] Fisher | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF], is about 150 µm for MIMAS MOX fuel. From [START_REF] White | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF] and [START_REF] Fisher | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF], the surface fraction of the Pu-rich spots is about 0.11 for MIMAS-ADU and 0.25 for MIMAS-AUC. A summary of these morphological information is presented in Table 1

Covariance function analysis for the spatial distribution of Pu-rich clusters

The covariance function is a tool that provides a clear and essential quantification of the spatial distribution of phases in a heterogeneous material (more details and other applications may be found in [START_REF] Jeulin | Random texture models for material structures[END_REF]; [START_REF] Lantuejoul | Geostatistical Simulations Models and Algorithms[END_REF]; [START_REF] Matheron | Random Set and Integral Geometry[END_REF][START_REF] Matheron | Eléments pour une théorie des milieux poreux[END_REF]; [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]; [START_REF] Torquato | Random Heterogeneous Materials[END_REF] and [START_REF] Wojtacki | Computing the elastic properties of sandstone submitted to progressive dissolution[END_REF]). Any binary image of the heterogeneous medium is one of the possible realizations generated by a specific stochastic process [START_REF] Lantuejoul | Geostatistical Simulations Models and Algorithms[END_REF][START_REF] Torquato | Random Heterogeneous Materials[END_REF], assumed to be stationary and ergodic. Each binary image, occupying the volume V , can be defined by the characteristic functions k X (x) of all the constitutive phases X

k X (x) = 1 if x ∈ X, 0 otherwise, (1)
With this definition, the volume fraction of phase X is the average of the characteristic function over the whole volume ϕ = k X (x) . The covariance function can be defined as the expectation of the autocorrelation of the characteristic function:

C(X, h) = E V k X (x)k X (x + h)dx = P {x ∈ X, x + h ∈ X}. (2)
In other words, the function C(X, h) gives the probability that both an arbitrary point of the phase X and its translation by the vector h belong to X. Thus, the covariance has the following properties:

C(X, h = 0) = |X|/|V | = ϕ and C(X, h → +∞) = (|X|/|V |) 2 = ϕ 2 , where h = h , | • | denotes the volume of a set.
Since the covariance function depends on the orientation of vector h , it is a convenient tool to investigate the anisotropy of a given microstructure. Six different 2D micrographs of MIMAS-ADU MOX microstructures have been extracted from the work of [START_REF] Largenton | MOX fuel effective behavior modeling by a micro-mechanical nonuniform transformation field analysis[END_REF]. These images have been binarized to extract the Pu-rich clusters and to generate the covariance maps of these inclusions. Figure 1 shows the averages over the six images of the covariances along three different directions (x, y and diagonal). Since the covariance is similar for all investigated directions, we assume in the sequel of the paper that the plutonium-rich clusters are isotropically distributed in the fuel. Moreover, the average information extracted from the covariance function C(X, h), where X is the Pu-rich where an overbar • denotes the average over the 6 different images.

Irradiation bubbles

Irradiation triggers the accumulation of rare gases such as Krypton or Xenon in Pu-rich clusters, which results in the apparition of irradiation bubbles (or pores). The size of these pores varies with their radial position in the pellet, increasing from the periphery towards the central region (Figure 2). The bubble size also depends on the local burn-up of the fuel: [START_REF] Noirot | Detailed characterisations of high burn-up structures in oxide fuels[END_REF] give histograms of bubble sizes in the center of the MIMAS MOX peripheral agglomerates for four distinct burn-ups. They show that the size of these bubbles increases with the burn-up (the peak of histograms increases from 2 µm to 13 µm with the burn-up). Moreover the largest bubbles are located in the central zone of the pellet, with a maximal size ranging from about 50 µm to 80 µm [START_REF] Noirot | Detailed characterisations of high burn-up structures in oxide fuels[END_REF][START_REF] Fisher | Measurement and analysis of fission gas release from BNFL's SBR MOX fuel[END_REF]. The pore volume fraction per cluster depends on their location and on the local burn-up rate as well (Figure 3): the porosity increases along the radius towards the central zone as well as with an increasing burn up rate. As can be seen in Figure 3, the maximal porosity found in clusters is 34%. A brief summary of this morphological information is given in Table 2 [START_REF] Noirot | Detailed characterisations of high burn-up structures in oxide fuels[END_REF].

Reconstruction of computational microstructures

The microstructure of MIMAS MOX is quite specific, strongly heterogeneous and evolutive with irradiation.

Simplifying assumptions concerning the geometry and the behavior of the constituent phases of this material are introduced hereafter in view of the computational modelling of representative volume elements.

Geometry

The geometry of irradiated MIMAS MOX has been idealized as a two-phase microstructure, with clusters of pores embedded in a surrounding matrix. The clusters, supposed to be spherical in shape and with a single size, were randomly distributed in the material. Spherical and single-sized pores were randomly distributed inside the clusters.

Taking into account the morphological information given in Section 2, the numerical samples used in the full-field simulations were reconstructed with the following set of parameters:

• porosity in the clusters: p c ∈ {0.1, 0.3}

• pore diameter: d p ∈ {4 µm, 8 µm, 16 µm, 32 µm} • volume fraction of clusters in the total volume: ϕ c = 0.15

• cluster diameter: d c = 100 µm
• number of clusters in the numerical samples2 : n c = 13.

Constitutive relations of the phases

For simplicity, by lack of experimental data for the individual phases, a single material behavior is attributed to the UO 2 -rich matrix (shown in black in Figure 4), to the coating phase in the case of MIMAS-ADU and to the Pu-rich spots (in grey in Figure 4). Following the aforementioned simplifications, a simple two-phase medium consisting of a homogeneous MOX matrix containing clusters of irradiation bubbles was investigated. It is generally admitted that the stationary creep of MOX fuel results from two creep mechanisms, namely diffusion creep and dislocation creep. The identification of material parameters for stationary creep of MOX fuel was done by [START_REF] Routbort | Compressive creep of mixed-oxide fuel pellets[END_REF]; [START_REF] Routbort | Creep of mixed-oxide fuel pellets at high stress[END_REF], and [START_REF] Roberts | Application of deformation maps to the study of in-reactor behavior of oxide fuels[END_REF]. The dependence of the creep rate on stress is described by the superposition of a linearly viscous flow and of a power-law creep relation. In the proposed models the maximal temperature in the range of validity is about T max ≈ 2000 K. A high temperature enhancement of the model was then proposed by [START_REF] Slagle | Fuel transient deformation[END_REF] who added a third term to increase the accuracy of the model up to T max = 2900 K. The final unidimensional constitutive relation of [START_REF] Slagle | Fuel transient deformation[END_REF] reads as follows:

εvp = A 1 d 2 e -Q1/(RT ) σ n1 + B 2 e -Q2/(RT ) + B 3 e -Q3/(RT ) σ n2 , (3)
where: εvp [1/s] is the strain rate, σ [MPa] is the stress, R [J/(mol • K)] denotes the gas constant (8.3144598

J/(mol • K)), d [µm]
is the grain size, the Q i [J/mol] denote activation energies, and A 1 , B 2 , B 3 depend on the density and oxygen-to-metal ratio. The parameters used in the present study are specified in Table 3. They correspond to a 0.97 density, a 1.97 oxygen-to-metal ratio and a 20 µm grain size (according to the first batch used in the initial study of [START_REF] Routbort | Compressive creep of mixed-oxide fuel pellets[END_REF]). According to [START_REF] Routbort | Compressive creep of mixed-oxide fuel pellets[END_REF], the first term in equation ( 3) with a linear dependence on the stress and an inverse-square-grain-size dependence is associated with a Nabarro-Herring stress-enhanced diffusion process which is likely to control the deformation in the linearly viscous flow region. In this regime, deformation is controlled by the diffusion of vacancies.

In the power law region, according again to [START_REF] Routbort | Compressive creep of mixed-oxide fuel pellets[END_REF], the deformation is due to the motion of dislocations (hence the nonlinear dependence on stress of the second term in equation ( 3)). A threedimensional version of Slagle's law has been implemented in CraFT (see Section 4). This law is used in the following to describe the viscoplastic behavior of the MOX matrix phase. The fission process also induces an athermal creep (whose mechanisms are not well known). Following [START_REF] Massih | Models for MOX fuel behavior[END_REF], this additional creep mechanism may be described by an additional term in equation (3) in the form C Ḟ σ, where C is a constant and Ḟ is the fission rate. Accounting for this term may be important to investigate fuel deformation under base irradiation conditions. This additional contribution to the strain rate can be incorporated into the first term of equation (3) (the linearly viscous term). The analytical model proposed in section 5 and based on the sum of two strain-rate potentials remains unchanged. In the present study focused on the behavior of the fuel during a postulated RIA, which lasts few seconds, this fission term is neglected. The pores (or bubbles) in the irradiated material are subjected to an internal pressure P b , due to fission gases accumulated in these cavities. A porous material with no internal pressure in the pores (P b = 0) is said to be drained, and, in the other case (P b = 0), the porous material is said to be saturated. The present study focuses on the overall viscoplastic behavior of a saturated porous material with an incompressible matrix. As demonstrated in [START_REF] Vincent | Porous materials with two populations of voids under internal pressure: I. instantaneous constitutive relations[END_REF], in the particular case of an incompressible matrix, the macroscopic behavior of the saturated material can be obtained from the macroscopic behavior of the drained material by a shift along the axis of hydrostatic stresses, Σ m (the macroscopic hydrostatic stress) being replaced by Σ m + P b . In other words when the voids are saturated, with pressure P b in the voids, and when the matrix is incompressible at the microscopic scale, it is sufficient to perform the scale transition for a drained porous material (P b =0) and then replace the overall hydrostatic stress Σ m by Σ m + P b . With this result in mind, only the drained case (P b = 0) is considered in the remainder of this study without loss of generality.

A 1 B 2 B 3 Q 1 Q 2 Q 3 d n 1 n 2 µm 2 • MPa -n1 • s -1 MPa -n2 • s -1 MPa -n2 • s -1 J/mol J/mol J/
As shown by fractographies (see [START_REF] Guérin | Microstructure evolution and in-reactor behaviour of MOX fuel[END_REF]), subdivision of grains occurs in the Pu-rich spots located at the pellet periphery and at mid-radius inside the fuel pellet. This subdivision of grains in very small grains is reminiscent of high-burnup structures observed in highly irradiated uranium dioxide fuels. By contrast, in the central area of the pellet, there is no significant change in grain size. As pointed out by [START_REF] Guérin | Microstructure evolution and in-reactor behaviour of MOX fuel[END_REF], the temperature is higher in the central area and the thermal diffusion is efficient enough to induce recovery of the accumulation of point defects, thus preventing subdivision of grains. In the present study, the material surrounding the cavities (without distinguishing between the UO 2 rich matrix, the Purich spots and, for MIMAS-ADU, the coating phase) is considered as isotropic. The isotropy assumption is legitimate when the cavities are larger than the surrounding grains. This is the case in the periphery of the pellet, after subdivision of the grains and when the cavities are larger than the grains in the Pu-rich spots. In other situations, for example in the central zone of the pellet, the cavities and the grains may have comparable sizes. Then the isotropy assumption is only motivated by the fact that the anisotropy of the crystals is moderate due to the cubic fluorite structure of PuO 2 and UO 2 .

Microstructure generation

Periodic numerical samples have been generated following the morphological properties defined in Section 3.1. The proposed algorithm proceeds in two steps (both steps generate monodisperse sets of objects):

• the Metropolis algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF] is used to arrange randomly n c non-overlapping spherical clusters that occupy a ϕ c fraction of the entire volume,

• then the same algorithm is used to distribute randomly inside each cluster n p non-overlapping spherical pores, that occupy a p c fraction of the cluster volume.

Finally, having a table with the positions and the radii of the clusters and of the pores, an image of the microstructure is generated at a desired spatial resolution. Figure 4 shows 2D examples of generated microstructures.

The periodic unit-cells used in the full-field simulations were generated using the above described methodology. A single initial spatial distribution of clusters has been used in all the generated microstructures, with The different parameters corresponding to the different microstructures are summarized in Table 4. The total volume fraction of the pores is f = p c ϕ c . Additionally, two simple microstructures were investigated in further simulations. These microstructures (labelled 1.0 and 2.0) are cubic unit-cells containing a single void with porosity f equal to 0.015 and 0.045 respectively. full-field simulations. The numerical method is briefly introduced in the next section, as well as the influence of the discretization (spatial resolution).

Micro. d c [µm] d p [µm] p c n p f 1.0 - - - 1 0.
4 Numerical method and influence of discretization

Integration of the constitutive equations

The constitutive equations of Slagle's model used in this study can be formulated for triaxial stress as

     ε = εe + εvp , σ = C : εe , εvp = 3 2 A1 d 2 e -Q 1 RT σ n1-1 eq σ d + (B 2 e -Q 2 RT + B 3 e -Q 3 RT )σ n2-1 eq σ d .
(4)

The strain ε has been decomposed into its elastic and visco-plastic parts ε e and ε vp , a dot over a variable denotes its time derivative, σ d is the deviatoric part of the stress σ, σ eq = 3 2 σ d : σ d is the equivalent von Mises stress and C denotes the elastic stiffness tensor.

The elastic moduli of the material are supposed to be isotropic, and thus fully described by the shear modulus µ and the bulk modulus k C = 2µK + 3kJ ,

where J and K are the usual fourth-order projectors on hydrostatic and deviatoric symmetric tensors of order 2. As a consequence the constitutive relations (4) can be further simplified as

     σd = 2µ εd - 3 2 A 1 d 2 e -Q 1 RT σ n1 eq σ d σ eq - 3 2 (B 2 e -Q 2 RT + B 3 e -Q 3 RT )σ n2 eq σ d σ eq , σm = 3k εm , (5)
where σ m = 1 3 Trσ denotes the hydrostatic component of the stress. Replacing the time derivatives ε and σ by the finite difference approximations

ε(t) ε(t) -ε(t -δt) δt and σ(t) σ(t) -σ(t -δt) δt (6)
and using a purely implicit scheme (backward Euler), the time-discretized constitutive relations (5) read as

     σ d -σ d (-δt) = 2µ ε d -ε d (-δt) -δt 3 2 A 1 d 2 e -Q 1 RT σ n1-1 eq σ d -δt 3 2 (B 2 e -Q 2 RT + B 3 e -Q 3 RT )σ n2-1 eq σ d , σ m -σ (-δt) m = 3k(ε m -ε (-δt) m ) , (7)
where ε d is the strain deviator (recall that the elastic part of the strain is compressible, while the viscous part is incompressible). For simplicity, the dependence in time is omitted for the current time, and thus ε and σ denote respectively the strain ε(t) and the stress σ(t) at time t, their values at time tδt being denoted by a superscript (-δt), ε (-δt) = ε(tδt) and σ (-δt) = σ(tδt).

Introducing a so-called "elastic trial" stress σ T defined by

σ T = σ (-δt) + C : ε δt
(this is the value that stress would have if no visco-plastic contribution occurred during tδt and t), one has

     σ d = σ T d -3µδt A 1 d 2 e -Q 1 RT σ n1-1 eq + (B 2 e -Q 2 RT + B 3 e -Q 3 RT )σ n2-1 eq σ d , σ m = σ T m . (8)
It follows from the first equation in (8) that σ d and σ T d are colinear. Thus, knowing the state of the stress and strain tensors at time tδt and knowing the strain at time t, σ T d can be evaluated straightforwardly and thus finding the stress at time t is reduced to a non linear equation whose unknown is σ eq

σ eq + 3µδt A 1 d 2 e -Q 1 RT σ n1 eq + 3µδt B 2 e -Q 2 RT + B 3 e -Q 3 RT σ n2 eq -σ T eq = 0 . (9)
In the present study, this scalar zero-finding equation has been solved using Muller's method [START_REF] Muller | A method for solving algebraic equations using an automatic computer[END_REF]. Once σ eq has been assessed, the whole tensor σ is obtained simply by:

   σ d = σ eq σ T eq σ T d σ m = σ T m . (10)
With this time-integration scheme, the constitutive equations ( 7) can be re-written symbolically with a function G σ

(x) = G ε(x), σ (-δt) (x), ε (-δt) (x) (11)
which makes the problem well suited to the use of the FFT-based homogenization method.

FFT-based homogenization method

The numerical approach that has to be chosen to carry out this study must take into account the very high complexity of the microstructures involved. Indeed, the dispersion of the clusters in the volume and the random distribution of porosities within the clusters, with the great disparity that exists between the size of the clusters and that of the pores, requires numerical simulations based on complex and large microstructural data. The computational method used here is the FFT-based homogenization method, introduced in Moulinec and Suquet (1994), which has long proven its capability to deal with this type of problem. The method, compared to the finite element method, offers the advantages of a lower computational cost and of being a meshless method -all the operations being performed on a regular grid. Its principles will be briefly recalled below, the interested reader is referred to the literature (among others Moulinec andSuquet (1994, 1998); [START_REF] Müller | Mathematical versus experimental stress analysis of inhomogeneities in solids[END_REF]; [START_REF] Eyre | A fast numerical scheme for computing the response of composites using grid refinement[END_REF] 2018)) for more information.

Lippmann-Schwinger equation

Let us consider a periodic unit-cell occupying a volume V of heterogeneous material submitted to an overall prescribed strain E. Periodicity conditions are imposed on the boundary ∂V of the unit cell. The properties of the phases are described by constitutive equations in the form σ = G ε, ... . Therefore, the local problem reads

         σ (x) = G ε(x), ... , div (σ (x)) = 0, ε (u (x)) = ε (u * (x)) + E, u * #, σ • n -#, (12)
where u * denotes the fluctuating part of the displacement. It satisfies periodicity conditions on the boundary of V (notation #), while the traction σ • n is antiperiodic (notation -#) for equilibrium purposes. By introducing a reference material with homogeneous stiffness c 0 , and the associated Green operator Γ 0 , the problem can be reduced to the Lippmann-Schwinger equation, which can be written in real space

ε (x) = -Γ 0 * τ (x) + E, (13)
where a polarization field τ has been introduced and defined as

τ (x) = σ(x) -c 0 : ε(x) = G ε(x), ... -c 0 : ε(x), (14)
and where * denotes a convolution operation. Relation (13) can be conveniently written in Fourier space as

ε (ξ) = -Γ 0 (ξ) :τ (ξ) ∀ξ = 0, ε (0) = E, (15)
where ξ denotes the angular frequency belonging to the reciprocal lattice associated with the unit-cell V .

Recall that a convolution in real space is replaced by a simple product in Fourier space.

Fixed-point algorithm

The fixed-point iteration method, proposed in [START_REF] Moulinec | A numerical method for computing the overall response of nonlinear composites with complex microstructure[END_REF] to solve the Lippmann-Schwinger equation ( 13), was adopted in the present study. The algorithm can be summarized by the following recursive relation:

ε k+1 (x) = -Γ 0 * σ k (x) -c 0 : ε k (x) + E σ k (x) = G ε k (x), ... ( 16
)
where k and k + 1 refer to two successive iterates. When the overall stress Σ is prescribed, the algorithm must be slightly modified into

ε k+1 = -Γ 0 * σ k (x) -c 0 : ε k (x) + E k+1 , σ k (x) = G ε k (x), ... , E k+1 = ε k -c 0 -1 : σ k -Σ (17)

Convergence tests

The algorithm requires a convergence test to decide when the iterating process must be stopped. The criterion proposed in [START_REF] Moulinec | A fast numerical method for computing the linear and the nonlinear mechanical properties of composites[END_REF] is based on an equilibrium test which consists in evaluating

err 1 (k) = div σ k 2 1/2 σ k (18)
and in comparing it, at each iteration, with the prescribed accuracy ζ 1 : the equilibrium test is satisfied when err 1 (k) < ζ 1 . In ( 18), the generic notation . has been used to denote the Frobenius norms of a vector v or a second-order tensor τ

v 2 = i=1,2,3 v 2 i , τ 2 = τ ij τ ij = i=1,2,3 j=1,2,3 τ 2 ij .
But as discussed in [START_REF] Moulinec | Convergence of iterative methods based on Neumann series for composite materials: Theory and practice[END_REF], this criterion imposes difficult conditions which is why the more relaxed criterion introduced by [START_REF] Monchiet | A polarization-based FFT iterative scheme for computing the effective properties of elastic composites with arbitrary contrast[END_REF] was used in this study. This criterion takes a simple form in Fourier space

err 1 (k) =          ξ 2 ξ. σ k (ξ) 2 ξ 2 - ξ. σ k (ξ).ξ 2 ξ 4 ξ σ k (ξ) 2          1 2 . ( 19
)
In addition, since the mechanical tests performed are creep tests or, more generally, tests where the macroscopic stress is imposed, a second convergence test must be performed at each iteration to determine if the macroscopic conditions prescribed are satisfied

err 2 (k) = σ k -Σ Σ < ζ 2 ( 20
)
where Σ is the prescribed overall stress and where ζ 2 is the requested accuracy for macroscopic conditions. In all the computations presented below, the convergence thresholds are ζ 1 = 10 -3 and ζ 2 = 10 -5 . In other words, iterations are stopped as soon as both conditions err 1 (k) < 10 -3 and err 2 (k) < 10 -5 are satisfied.

Stabilization of the overall strain rate

The mechanical conditions of a creep test are simulated. The macroscopic stress is prescribed and the response of the material is described by the evolution of the overall deformation with time. After a transient regime whose duration varies with temperature and mechanical loading and, to a lesser extent, with the sample considered, the overall strain rate stabilizes and remains constant.

In the present study, the overall strain rate is estimated by a five point stencil

Ė(t) -E(t + 2δt) + 8E(t + δt) -8E(t -δt) + E(t -2δt) 12 δt . ( 21
)
The stabilized regime is reached -and the loading process is then stopped -when the relative evolution of the strain rate is small. The stopping criterion used is

Ė(t) -Ė(t -Kδt) < η L Ė(t) . (22)
In the present work, the typical value for K is 10 (which means that two strain rates are compared with 10 loading steps of difference), and the required accuracy is η L = 10 -5 .

Spatial discretization

The question of spatial discretization is of crucial importance in determining the ability of numerical simulations to capture local variations in the mechanical fields involved and hence to accurately describe the effective properties.

The FFT-based method imposes a spatial discretization along a regular grid, the basic element of which is a voxel. The discretization study should proceed by performing several simulations of the same structure at different resolutions (i.e. by varying the numbers of voxels in the volume considered). Its cost would be beyond our currently available computational means. The approach adopted here was to study a simplified microstructure that is nevertheless representative of the actual geometry. This simplified microstructure, shown in Figure 7, consists of a single cluster located at the center of the image, with a volume fraction of 15%, in which 10 pores have been placed randomly. The two porosities of interest, 10% and 30% of pores per cluster, are investigated. The mechanical loading is a uniaxial creep test in compression (Σ = -Σ e 1 ⊗ e 1 ), with the following values of the thermomechanical parameters T = 2473 K and Σ=100 MPa. f = 1.5% Ė11 (10 -6 1/s) 6.16 6.14 6.13 6.11 6.11 err rel [%] 0.9 0.6 0.3 0.05 f = 4.5% Ė11 (10 -6 1/s) 7.41 7.30 7.23 7.18 7.16 err rel [%] 3.54 2.01 0.94 0.3 -Table 5: Centered cluster with 10 pores. Overall strain rates and relative differences (with respect to the finest discretization) as a function of the discretization along the pore diameter.

Several discretizations of the unit-cell were examined, 32 3 , 64 3 , 128 3 , 256 3 and 512 3 voxels, corresponding to an increasingly refined resolution along the diameter of a single pore. Table 5 summarizes the results obtained. The pore volume fraction of 4.5% (30% of pores in a cluster) is more sensitive to discretization than the porosity 1.5%. A compromise between accuracy and computational cost was achieved by selecting a resolution of 29 and 42 voxels per diameter of pore for 1.5% and 4.5% porosity respectively. The precision on the overall strain rate is then less than 0.1% and 0.5% respectively.

Full-field simulations

In order to model the effect of pores on the overall behaviour of the irradiated MOX fuel, 10 different temperatures and 10 stress levels were investigated in the simulations:

Temperature: T [K] ∈ {1573, 1673, ..., 2473} , Stress: Σ [MPa] ∈ {10, 20, ..., 100} .
The simulations were mostly focused on the uniaxial compressive creep test, given by:

Uniaxial: Σ = Σ (-e 1 ⊗ e 1 ) . (23)
Due to the chosen discretization and related computational cost, the simulations performed on the most refined microstructures were limited (see Table 6). Additional creep tests, under hydrostatic and axisymmetric deviatoric loadings, were carried out in order to adjust the semi-analytical model presented in Section 5. These additional simulations were performed at temperature T = 1973 K and at only one stress level for all generated microstructures. The hydrostatic and deviatoric loadings are given respectively as:

Hydrostatic: Σ = Σ (-e 1 ⊗ e 1 -e 2 ⊗ e 2 -e 3 ⊗ e 3 ) , Deviatoric: Σ = Σ -2 3 e 1 ⊗ e 1 + 1 3 e 2 ⊗ e 2 + 1 3 e 3 ⊗ e 3 . ( 24
)
Detailed information about the microstructures and the different loading conditions are summarized in Table 6.

It is worth mentioning that the computations on the largest microstructure (micro. 1.4) were run with approximatively 4TB of RAM memory and 1296 cores.

Analytical model

Strain-rate potentials

An estimate of the effective viscoplastic behavior of the porous MIMAS MOX fuel can be proposed using a micromechanical approach. It is noted that the matrix material is governed by a strain-rate potential exhibiting two temperature-activated creep mechanisms Table 6: Summary of chosen geometrical and computational information related to all microstructures used in full-field simulations: overall porosity, number of pores per cluster, RVE edge length, number of performed simulations for uniaxial, hydrostatic and deviatoric loadings respectively. The second sub-column of N c uni is the number of investigated temperatures and stress levels.

ε = ∂ψ ∂σ (σ), ψ(σ) = ψ 1 (σ) + ψ n (σ), ( 25
) Micro. f [%] n p L [vox] N c uni N c hyd N c dev 1.0 1.
with ψ 1 (σ) = 1 2 Aσ 2 eq , ψ n (σ) = 1 n + 1 Bσ n+1 eq , (26)
where RT ) . Similarly the overall strain rate of the porous material also derives from an effective strain-rate potential (Ponte Castañeda and Suquet, 1997)

A = A 1 d 2 e -Q1/(RT ) , B = B 2 e -Q2/(RT ) + B 3 e -Q3/(
Ė = ∂Ψ ∂Σ (Σ) with Ψ(Σ) = Inf σ ∈ S (Σ) ψ(σ) , (27)
S (Σ) = {σ ∈ H , div σ(x) = 0 in V, σ.n -#, σ = Σ}. ( 28
)
Define two partial effective potentials as

Ψ 1 (Σ) = Inf σ ∈ S (Σ) ψ 1 (σ) and Ψ n (Σ) = Inf σ ∈ S (Σ), ψ n (σ) , (29)
and let σ, σ 1 , σ n be the local stress fields solution of the variational problems ( 27) and ( 29). Then, by the variational properties ( 29)

ψ 1 (σ) ≥ ψ 1 (σ 1 ) , ψ n (σ) ≥ ψ n (σ n ) ,
and therefore

Ψ(Σ) = ψ 1 (σ) + ψ n (σ) ≥ ψ 1 (σ 1 ) + ψ n (σ n ) = Ψ 1 (Σ) + Ψ n (Σ)
In other words Ψ 1 + Ψ n is a lower bound for the effective potential Ψ

Ψ(Σ) ≥ (Ψ 1 + Ψ n )(Σ). ( 30
)
Remark: Homogenization is a non-additive procedure, in the sense that the effective potential corresponding to the sum of two potentials (here ψ 1 and ψ n ) is not the sum of the individual effective potentials (here Ψ 1 and Ψ n ). This is a well-known fact in viscoelastic composites which gives rise to long memory effects, even when the individual constituents exhibit only short-term memory effects. These complex interactions between the two potentials will be neglected in the sequel, but the inequality (30) gives a rigorous status to this approximation. Note that [START_REF] Monerie | Overall viscoplastic behavior of non-irradiated porous nuclear ceramics[END_REF] approached the problem differently (and heuristically) by interpolating Ψ between Ψ 1 and Ψ n with a weight factor depending on the overall stress Σ.

Bounding the overall strain rate

Although a lower bound on the potential does not always result in a lower bound for its derivatives, in the present case the lower bound (30) implies a lower bound for one of the components of the overall strain rate

Ė = ∂Ψ ∂Σ (Σ). (31)
Recall that σ, σ 1 , σ n be the local stress fields solution of the variational problems ( 27) and ( 29). Euler's theorem for positively homogeneous functions, applied to ψ 1 and ψ n respectively, yields

       ∂Ψ 1 ∂Σ (Σ) : Σ = 2Ψ 1 (Σ) = 2 ψ 1 (σ n ) ≤ 2 ψ 1 (σ) = ∂ψ 1 ∂σ (σ) : σ , ∂Ψ n ∂Σ (Σ) : Σ = (n + 1)Ψ n (Σ) = (n + 1) ψ n (σ n ) ≤ (n + 1) ψ n (σ) = ∂ψ n ∂σ (σ) : σ . (32)
Adding up the two inequalities in (32), one gets that

∂Ψ 1 ∂Σ (Σ) : Σ + ∂Ψ n ∂Σ (Σ) : Σ ≤ ∂ψ 1 ∂σ (σ) + ∂ψ n ∂σ (σ) : σ = ε : σ , (33)
and finally using Hill's lemma Ė :

Σ ≥ ∂Ψ 1 ∂Σ (Σ) + ∂Ψ n ∂Σ (Σ) : Σ. (34)
Therefore, replacing the exact effective potential Ψ with its lower bound Ψ 1 + Ψ n yields a rigorous lower bound for the overall strain rate Ė in the direction of the applied stress Σ.

5.3 Estimate using the model of Leblond-Perrin-Suquet (LPS) [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF] It remains to bound or estimate the two partial effective potentials Ψ 1 and Ψ n . Following [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF], using the fact that ψ 1 and ψ n are power-law potentials, Ψ 1 and Ψ n can be written as

Ψ 1 (Σ) = ψ 1 (Λ 1 ) , Ψ n (Σ) = ψ n (Λ n ) , (35)
where Λ 1 and Λ n are positively homogeneous functions of degree 1 with respect to Σ. To simplify notations, ψ 1 and ψ n , which were introduced in (26) as functions of the tensorial variable σ, are considered in (35) as power-law functions of the single scalar variable σ eq . Λ n (Σ) is obtained by solving the equation [START_REF] Leblond | Exact results and approximate models for porous viscoplastic solids[END_REF])

F(Q, M ) = 0, (36)
where

F(Q, M ) = (1 + 2 3 f )Q 2 + f h(M ) + n + 1 n -1 1 h(M ) -1 - n -1 n + 1 f 2 , (37)
Q = Σ eq Λ n , M = Σ m Λ n , h(M ) = 1 + 1 n 3 |M | 2 n+1 n n , (38)
Λ 1 is found from eq. ( 36) with n = 1 and f denotes porosity,

Λ 1 = 9 4 f 1 -f Σ 2 m + 1 + 2 3 f 1 -f Σ 2 eq 1/2 . ( 39
)
There is no expression in closed-form for Λ n when n = 1 and the equation (36) has to be solved numerically. This can be done by a Newton-Raphson method, initialized by the approximation provided by the hollow sphere model of [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF]

(see Appendix A).

The strain rate predicted by the effective potential Ψ = Ψ 1 + Ψ n reads as

Ė = ψ 1 (Λ 1 ) ∂Λ 1 ∂Σ + ψ n (Λ n ) ∂Λ n ∂Σ . ( 40
)
∂Λn ∂Σ is obtained by taking the derivative of (36) with respect to Σ and the resulting expression is

∂Λ n ∂Σ = Λ n ∂F ∂M Σ m + ∂F ∂Q Σ eq 1 3 ∂F ∂M i + 3 2 ∂F ∂Q Σ d Σ eq . (41)
The lower bound (34) can be given a simple form. Using Euler's relation for the positively homogeneous functions Λ 1 and Λ n , one gets

∂Λ 1 ∂Σ : Σ = Λ 1 , ∂Λ n ∂Σ : Σ = Λ n , (42)
and therefore Ė :

Σ = ψ 1 (Λ 1 ) Λ 1 + ψ n (Λ n ) Λ n = AΛ 2 1 + BΛ n+1 n . ( 43
)
6 Full-field simulations versus analytical model

Modified porosity

In order for the comparison between full-field simulations and the analytical models developed in section 5 to be meaningful, a consistent definition of porosity has to be adopted on both sides. This is the role played in the GTN model by the parameters q 1 , q 2 , q 3 introduced by [START_REF] Tvergaard | On localization in ductile materials containing spherical voids[END_REF]. Roughly speaking these parameters account for the fact that in the model, the porous material will loose its carrying capacity when f = 1, whereas in the numerical simulation, the pores will percolate much earlier. Therefore the results of the full-field simulations at porosity f should be compared to the analytical model with a modified porosity f * = qf where q is a parameter between 1 and 2 (when the voids are spherical a classical value for the parameter q is 1.25). A similar notion of modified porosity has been used in several different studies (see [START_REF] Fritzen | Computational homogenization of elasticplastic porous materials[END_REF]Vincent et al., 2014, among others).

In the present study, the q parameter was identified from two simulations performed at T = 1973 K for all available microstructures (at each void volume fraction) under hydrostatic and axisymmetric deviatoric loadings in the form (24), with Σ = 30 MPa (hydrostatic loading) and Σ = 45 MPa (axisymmetric deviatoric loading).

The q parameter in the model was adjusted from its usual value for spherical voids (q = 1.25) to give a reasonable agreement between the model and the simulations (different microstructures at each void volume fraction) for both loadings and the following values were found q = 1.28 when f = 0.015, q = 1.35 when f = 0.045.

Then these values of q were used, with no further fitting, to simulate uniaxial compression tests at different stress levels and different temperatures. Comparisons with full-field simulations are shown in Figures 10 and11. The agreement is seen to be quite satisfactory over the whole range of stress and temperature investigated.

Discussion

From the results shown in Figures 10 and11, the following observations can be made. Figure 9: Identification of the parameter q using deviatoric and hydrostatic loadings. Comparison between the model (40) (solid line) and full-field simulations (full circles corresponding to the 3 different microstructures: 2.1, 2.2 and 2.3). Void volume fraction f = 0.045. Ėmatrix eq denotes the equivalent strain rate in the fully dense material.

• As expected, the ratio Ė11 / Ėmatrix 11 is always above 1, which is a clear indication that the cavities enhance the overall strain rate and soften the material.

• The effect of cavities increases with the stress: for example, when f = 0.045 and T = 2473K, the model predicts that the ratio Ė11 / Ėmatrix 11 increases from 1.21 to 1.53 when the stress increases from 10 to 100 MPa. It shows that in the high-stress regime (mainly associated with the non-linear term in the viscoplastic law of the matrix) the material is more sensitive to cavities than in the low-stress regime (mainly associated with the linear term in the viscoplastic law of the matrix).

• The evolution of the ratio Ė11 / Ėmatrix 11 with respect to the stress is highly nonlinear: a strong increase of this ratio is observed in a transition zone between a lower and an upper plateau where this ratio seems to be stabilizing. The effect of cavities saturates at high stress or low stress. The slope in this transition zone increases with temperature. Since the temperature strongly enhances the nonlinear term in the viscoplastic behavior, this observation is in line with the previous comment: an increase in temperature implies an increase in the nonlinear term, and thus, the transition zone from the low stress to the high stress regime is shortened.

• The effect of porosity on the ratio Ė11 / Ėmatrix 11 is quite significant. For example, when Σ = 100 MPa and T = 2473 K, both the simulations and the model show that this ratio increases from 1.14 to 1.53 when the porosity goes from 1.5 % to 4.5 %. This increase is less pronounced in the low stress-low temperature regime: when Σ = 10 MPa and T = 1573 K, the model predicts an increase in this ratio from 1.03 to 1.12 when the porosity goes from 1.5 % to 4.5 %. Consequently the high stress-high temperature regime (mainly associated with the nonlinear term in the viscoplastic behavior) is more sensitive to the cavities than the low stress-low temperature regime (mainly associated with the linear term in the viscoplastic behavior).

• The influence of porosity on the overall strain rate is significant: in the range of porosities, temperatures and stresses studied here, it can induce a strain rate up to 50 % higher than in the dense case (i.e. without porosity). We thus conclude that the presence of cavities plays a key role in the viscoplastic behavior of this material. In other words, taking into account the porosity due to irradiation when studying the viscoplasticity of this irradiated MOX fuel is of prime importance.

• In the case of a 1.5 % porosity, the changes in the strain rates with respect to the microstructures (1.1, 1.2, 1.3) are quite limited: as an example, when Σ = 100 MPa and T = 2473 K, the relative difference in strain rate between the microstructures 1.1 and 1.3 is only 0.3 %. It is more pronounced but still weak in the 4.5 % porosity case: when Σ = 100 MPa and T = 2473 K, the relative difference in strain rate between the microstructures 2.1 and 2.3 is 2.3 %. It may be concluded that moving from a microstructure with a small number of pores inside each cluster to a microstructure with a huge number of pores inside each cluster without modifying the porosity (it induces a decrease in the size of each cavity) has only a limited impact on the overall strain rate. This conclusion is strongly linked to the specific modeling of the material studied here under the physical conditions specified above. Note that, since the size of the cavities considered here are above a few microns, no size effect has been considered in the modeling: as a general rule, surface effects are predominant for nanosized inclusions and vanish when their size increases. Moreover, differences of about 5 % can be observed in the strain rate between the microstructure 2.0 (relative to voids distributed according to a cubic lattice in the matrix) and the other microstructures. From these comments, it may be concluded that the void clustering has only a rather limited impact on the overall strain rate.

Conclusion

MIMAS MOX fuels exhibit a very specific microstructure, strongly heterogeneous and evolutive with irradiation. Fresh fuels contain Pu-rich clusters. During irradiation, accumulation of rare gases generated by the fission triggers the onset of irradiation bubbles within the material. The Pu-rich clusters become highly porous, with pressurized bubbles. Concerning the Pu-rich clusters and their irradiation bubbles, data on their sizes and surface fractions have been collected from a literature review. A covariogram analysis for the Pu-rich clusters has been performed based on available digital images of fuels, showing a good agreement with the collected data.

Then, the geometry of this material has been idealized as a two-phase microstructure, a matrix in which spherical pores are gathered into spherical clusters. The matrix is governed by a strain-rate potential exhibiting two creep mechanisms and the pores are subjected to a homogeneous internal pressure (irradiation bubbles). Only the secondary creep regime was addressed here.

Periodic cells have been generated with a single spatial distribution of 13 clusters in which pores are randomly distributed. Two porosities have been considered together with four different pore diameters. The microstructure with the smallest pore size contains more than 1500 pores in each cluster.

These periodic cells have been discretized into a regular cubic array of voxels and FFT-based numerical full-field simulations have been performed. First of all, the nonlinear behavior of the matrix has been integrated in time by a fully implicit algorithm which has been implemented in the FFT software. Then, a discretization study has been performed on a simplified cell with a single porous cluster. The spatial resolution for the different periodic cells required to achieve mesh independence ranges from 324 3 to 2592 3 voxels. In terms of degrees of freedom, the problem to be solved for the highest resolution amounts to several tens of billions of dof's. These massive full-field simulations have been performed on a super-computer, with up to 1296 cores and 4 TB of RAM.

An analytical estimate for the effective behavior of this porous viscoplastic material has been derived through a variational approach. A modified porosity has been introduced in the model and identified from two simulations for all microstructures at a given temperature under hydrostatic and deviatoric loadings. Then, the model has been used to simulate the uniaxial compression tests at different stress levels and different temperatures. Comparisons with full-field simulations show a good agreement over the whole range of stress and temperature investigated.

As expected, it has been observed that the presence of cavities leads to an increase in the overall strain rate and thus softens the material. Moreover, it seems that the nonlinear term in the viscoplastic law is more sensitive to cavities than the linear term. The effect of the cavities tends to saturate at high stress or low stress and the transition zone between these two regimes depends on the temperature. From the full-field simulations, it may be concluded that gathering the voids into clusters has only a rather limited impact on the overall strain rate of the voided material. But the global influence of the cavities on the overall strain rate is significant: a strain rate up to 50 % higher than in the dense case (i.e. without porosity) has been obtained. It can be concluded that porosity plays a key role in the creep behavior of this material. f 2 = q 2 f are introduced. These modified porosities enter the prefactors of the mean overall stress and the equivalent overall stress respectively:

             Σ 1 = 9 4 f 1 Σ 2 m + (1 + 2 3 f 2 )Σ 2 eq 1/2 , Σ n =   9 4 1 -f 1 n(f -1 n 1 -1)

  covariance of Pu-rich clusters

Figure 1 :

 1 Figure 1: Distribution of Pu-rich clusters in MIMAS-ADU MOX : image and covariance function. The Pu-rich clusters are in red.

Figure 2 :

 2 Figure 2: Apparition of irradiation bubbles in the Pu-rich clusters. The bubble size changes with distance from the pellet center (MIMAS MOX Noirot et al. (2008)).

Figure 3 :

 3 Figure3: Dependence of porosity in the Pu-rich agglomerates on the local burn-up for MOX fuel[START_REF] Noirot | Detailed characterisations of high burn-up structures in oxide fuels[END_REF].

  (a) nc = 5, ϕc = 0.45, np = 50, pc = 0.2 (b) nc = 10, ϕc = 0.5, np = 100, pc = 0.25

Figure 4 :

 4 Figure 4: Example of periodic unit-cells: pores (white), clusters (grey) and matrix (black).

Figure 5 :

 5 Figure 5: Pu-rich clusters distribution used to generate the periodic cells for the full-field simulations.

Figure 6 Figure 6 :

 66 Figure 6 shows one cluster of each microstructure. All seven generated samples are used in the further

  ; Zeman et al. (2010); Kabel et al. (2014); Willot et al. (2014); Schneider et al. (2016); Moulinec et al. (

  (a) pore volume fraction f = 0.015 (b) pore volume fraction f = 0.045

Figure 7 :

 7 Figure 7: Pores distribution in the unit-cell used in the study of discretization.

Figure 8 :

 8 Figure8: Identification of the parameter q using deviatoric and hydrostatic loadings. Comparison between the model (40) (solid line) and full-field simulations (full circles corresponding to the 4 different microstructures: 1.1, 1.2, 1.3 and 1.4). Void volume f = 0.015. Ėmatrix eq denotes the equivalent strain rate in the fully dense material.

Figure 10 :Figure 11 :

 1011 Figure 10: Uniaxial compression. Void volume fraction f = 0.015. between the model (40) (solid line) and full-field simulations (full circles) performed with the FFT-based method. Microstructure 1.0 (black circle), 1.1 (red), 1.2 (green), 1.3 (blue), 1.4 (grey). Ėmatrix 11 denotes the equivalent strain rate in the fully dense material.

  Figure A.12: Deviatoric hydrostatic loadings. Comparison between the model (A.16) (solid line) and full-field simulations (full circles corresponding to the 4 different microstructures: 1.1, 1.2, 1.3 and 1.4). Void volume fraction f = 0.015. Ėmatrix eq

  Figure A.14: Uniaxial compression. Void volume fraction f = 0.015. Comparison between the model (A.16) (solid line) and full-field simulations (full circles) performed with the FFT-based method. Microstructure 1.0 (black circle), 1.1 (red), 1.2 (green), 1.3 (blue), 1.4 (grey). Ėmatrix 11 denotes the equivalent strain-rate in the fully dense material.

  Figure A.15: Uniaxial compression. Void volume fraction f = 0.045. Comparison between the model (A.16) (solid line) and full-field simulations (full circles) performed with the FFT-based method. Microstructure 2.0 (black circle), 2.1 (red), 2.2 (green), 2.3 (blue). Ėmatrix 11 denotes the equivalent strain-rate in the fully dense material.

Table 1 :

 1 .

	Descriptor	Value
	Surface fraction 0.11 -0.25
	Maximal size	150 µm

Summary of morphological information concerning the Pu-rich clusters (agglomerates) in MIMAS MOX fuel.

Table 2 :

 2 . Summary of morphological information concerning the irradiation bubbles.

	Descriptor	Value
	porosity in the clusters p c	0 -0.34
	characteristic size	2 µm -80 µm

Table 3 :

 3 Parameters of Slagle's law (3) used in the present study.

	mol	µm	-	-

Table 4 :

 4 Summary of geometrical information of the microstructures used in the full-field simulations.

	015

http://craft.lma.cnrs-mrs.fr/

The number of clusters is limited by the available RAM memory.

Acknowledgments

Some of the full-field simulations presented here were made possible thanks to the resources of the CCRT (Centre de Calcul Recherche et Technologie) of the CEA/DIF. HM and PS have received funding from Excellence Initiative of Aix-Marseille University -A*MIDEX, a French "Investissements d'Avenir" programme in the framework of the Labex MEC.

Appendix A The hollow sphere model of [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF] Michel and [START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF] have noticed that the two effective potential Ψ 1 and Ψ n can be written as

where Σ 1 and Σ n are positively homogeneous functions of degree 1 of Σ such that [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF] Σ n ≥ Σ p ∀n ≥ p.

(A.

2)

The effective stresses Σ 1 and Σ n are related to Λ 1 and Λ n through

Combining ( 30) and (A.2) yields a first simple and rigorous lower bound for Ψ:

It can easily be checked that the lower bound (A.4) is nothing else than the variational lower bound obtained by means of an isotropic linear elastic comparison material [START_REF] Ponte Castañeda | The effective mechanical properties of nonlinear isotropic composites[END_REF][START_REF] Suquet | On bounds for the overall potential of power law materials containing voids with an arbitrary shape[END_REF][START_REF] Ponte Castañeda | Nonlinear composites[END_REF]. Σ 1 is related to the linear effective properties of a linear porous material with the same microstructure as the original nonlinear one, but with an incompressible matrix with a shear modulus µ 0 through

When the pores are randomly distributed (shape, orientation and location are random), the Hashin-Strikman upper bound can be used to bound M and therefore Σ 1 from below

and a rigorous lower bound for Ψ can be obtained by combining (A.4) with (A.6). However, as is well-known from other studies, the variational lower bound (A.4), where Σ n is replaced by its lower bound Σ 1 , is rather sharp when the overall stress is deviatoric, but rather poor under hydrostatic loading. We therefore use the sharper lower bound

with a more accurate estimate for Σ n . [START_REF] Michel | The constitutive law of nonlinear viscous and porous materials[END_REF] have proposed to use the exact response of a hollow sphere under hydrostatic loading to get an estimate for Σ n which is more accurate under hydrostatic stress :

The strain-rate predicted by this sharper lower bound reads as

where ∂Σ1 ∂Σ and ∂Σn ∂Σ can be deduced from (A.6) and (A.9) respectively

where it is recalled that Σ d denotes the overall stress deviator. Using Euler's relation for the positively homogeneous functions and Σ n , one gets

and therefore

For uniaxial compression with Σ = -Σe 1 ⊗ e 1 (Σ ≥ 0) one gets

A first comparison between full-field simulations and the hollow sphere with the actual porosity shows that the model with the actual porosity f underestimates the overall strain-rate. Conversely, introducing the modified porosity f * = qf and applying the same methodology of identification of the q parameter as in section 6.1, leads to an overestimation of the overall strain rate for uniaxial loading over the whole range of temperature. Therefore, two parameters q 1 and q 2 , corresponding to two modified porosities f 1 = q 1 f and The corresponding effective relation under triaxial stress reads as

Under hydrostatic loading Σ = -Σ(e 1 ⊗ e 1 + e 2 ⊗ e 2 + e 3 ⊗ e 3 ) (Σ ≥ 0) one gets

Σ which leads to the following estimation of Tr Ė:

(A.17)

Similarly, under deviatoric loading with Σ = Σ eq (-2 3 e 1 ⊗ e 1 + 1 3 e 2 ⊗ e 2 + 1 3 e 3 ⊗ e 3 )

Σ eq which gives:

Finally, for uniaxial compression with: Σ = -Σe 1 ⊗ e 1 (Σ ≥ 0) one gets:

Σ and the overall strain-rate in the direction of the compressive stress reads as:

The parameters q i are identified for each microstructure separately. The parameter q 1 is found from the hydrostatic creep test, the parameter q 2 from the uniaxial test at T = 2473 K and under the largest applied stress Σ = 100 MPa. Finally the obtained parameters are averaged among different microstructures, at each value of porosity.

The predictions of the calibrated model (by identification of q i ) are compared with full-field simulations in Figures A.14 and A.15. The agreement is seen to be quite satisfactory.