Homogenization approach to model the thermal-mechanical behavior of melting fuel material

Abstract For nuclear pressurized water reactors (PWR), safety analysis are carried out on operating modes involving quick and large amplitude changes in thermal power produced. In this frame, assessments have to conservatively account for initiating fuel melting localized in the pellet central region, where temperatures and stresses are the highest. In the present paper we propose a thermal-mechanical coupled model to describe the behavior of fuel pellet at melting onset, developed by means of homogenization methods. The model has been implemented in the finite element code Cast3m and a simple test case is presented to verify and validate the modeling principle. The dependency of the thermal mechanical behavior of the system on the melted fraction is discussed.

For nuclear pressurized water reactors (PWR), safety analysis are carried out on operating modes involving quick and large amplitude changes in thermal power produced. In this frame, assessments have to conservatively account for initiating fuel melting localized in the pellet central region, where temperatures and stresses are the highest. In the present paper we propose a thermal-mechanical coupled model to describe the behavior of fuel pellet at melting onset, developed by means of homogenization methods. The model has been implemented in the finite element code Cast3m and a simple test case is presented to verify and validate the modeling principle. The dependency of the thermal mechanical behavior of the system on the melted fraction is discussed.

Introduction: power transients and the onset of fuel melting
The operation of French nuclear power plants is nowadays evolving towards a more flexible functioning to follow grid demand. Flexibility of power production implies that nuclear fuel rods have to face quick and large amplitude changes in the thermal power produced by the plant. Those power transients induce strong thermal-mechanical stresses in the fuel element. Safety analyses are carried out for those kind of operations to define technical specification that guarantee the integrity of the cladding.
In this study we consider a hypothetical power transient leading to a maximum linear power produced in the fuel much higher than normal operating conditions. This power leads to the onset of solidliquid melting transition where temperature locally overcomes the melting threshold. To study the macroscopic behavior of the melting fuel material, we develop a model based on recent homogenization techniques that describes its thermal-mechanical behavior. This article focuses on the validation of the modeling principle before its implementation in ALCYONE, a multi-dimensional fuel performance code co-developed in the PLEIADES platform by CEA, EDF and FRAMATOME [1]. ALCYONE is a complex multi-physics numerical tool, used to simulate the fuel element behavior during normal and off-normal operation conditions. This requires the coupling of the neutronic, thermohydraulic, thermo-mechanical and physicochemical models. A simplified flowchart of ALCYONE is presented in A. To validate the proposed melting model, we implement it in the finite element code Cast3m [2], verify its convergency and test its calculated outputs in a simple case. All the aspects related to material knowledge are treated by means of the code generator MFRONT, which is developed within the PLEIADES platform [3,4]. Its future implementation in ALCYONE is meant to evaluate the effect on clad deformation induced by the progressive melting of the fuel and to guide the development of a measurement method adapted to study the melting scenario in a dedicated instrumented device in the JHR material testing reactor [5].
In this paper, we first present in section 2 a review of the thermal-mechanical behavior of solid fuel material. For this specific work, we use the model proposed by Monerie and Gatt in Ref. [6]. Section 3 is dedicated to the thermal-physical aspects of the melting process. In section 4 we proceed to the description of the analytic model proposed, here stating the methodology and hypothesis used and detailing the effects on the mechanical behavior.
In section 5 we discuss the resulting behavior obtained by the implementation of the model in Cast3M. Conclusions and further work are presented in section 6.

The thermal-mechanical behavior of the solid fuel material
Dense uranium dioxide presents at high temperature an elastoviscoplastic behavior. Frost and Ashby in Ref. [7], describe that uranium dioxide presents two predominant creep mechanisms which are activated above 1000 C and are strongly temperature dependent: the deformation map presented in Fig. 1 shows that the oxide undergoes dislocation diffusion mechanism (power law creep in the figure) when submitted to high stresses/high temperatures, it instead goes through cation diffusion when lower stresses/lower temperatures are applied.
Dislocation-creep and cation diffusion-creep mechanisms for an incompressible material can be described by means of the Norton potential j vp i ðs; n i Þ: where s eq is the second invariant of the stress (or Von Mises stress) s eq ¼ ffiffiffiffiffiffiffiffiffiffiffi ffi 3 2 s : s q , s being the deviatoric component of the stress tensor s, n i represents the Norton exponent for the i th mechanism, dislocation or cation diffusion, s 0 and _ ε 0 represent material constants. The strain rate for the specific mechanisms is consequently defined as the derivative with respect to the stress tensor of the total dissipative potential j vp i ðs; n i Þ [8]: The solid fuel material presents by fabrication a porous fraction, which evolves during irradiation as a consequence of fission phenomena and impacts the fuel material behavior. The porous fraction present by fabrication is meant to decrease during irradiation due to the densification process taking place during the first irradiation phases. A second type of porosity appears during its life in the reactor, due to gaseous products resulting from nuclear fission events. The presence of porosities leads to the appearance of a compressible deformation in the viscoplastic nuclear fuel. Several studies have been carried out in the past to define the mechanical behavior of viscoplastic matrices containing voids. Rice and Tracey [9] first evaluated the deformation of a stiff-plastic matrix presenting a single void, whose result was then introduced in the Gurson's [10] derivation of the yield criterion of rigid-plastic materials containing arbitrary void fraction. Interaction effects were then introduced by Tvergaard [11] according to numerical simulations and successively revised by Perrin and Leblond in 1994 [12].
For non linear viscous matrices, as the nuclear fuel, the studies on infinite non linear matrix presenting a single void were first carried out by Budiansky et al., in 1982 [13] and then addressed by Duva andHutchinson in 1984. Duva in 1986 [14] examined the numerical problem for an arbitrary number of voids and proposed to use Gurson's potential as a basis for the viscous strain rate potential. Bounds for the potential were derived in 1988 by Ponte Castaneda and Willis [15] by means of variational principle. Finally in 1992, Michel and Suquet [16] derived a simple model to find the exact solution of a hollow sphere under hydrostatic pressure, for incompressible non linear behavior as well as porous one. On this basis, the thermal-mechanical modeling proposed by Monerie and Gatt [6] describes the behavior of a non-linear compressible viscoplastic medium which porous fraction, f, is randomly distributed and the material can be considered isotropic. This is expressed as a function of the creep mechanisms n i and of the porous fraction f thorough the dissipative elliptical potential j vp i ðs; f ; n i Þ: where s 0 and _ ε 0 are the reference stress and strain rate respectively; s H ¼ 1 3 trðsÞ is the hydrostatic part of stress tensor. An exhaustive discussion about the choice of parameters Aðf Þ and Bðf Þ for UO 2 can be found in Ref. [6]. In the application cases proposed in section 5 we retain the form deduced by the selfconsistent estimate by Michel and Suquet [16].
Monerie in Ref. [6] introduces a temperature dependent function qðTÞ to couple the two creep mechanisms. It represents the importance of each mechanism depending on temperature: it equals 0 when only diffusion creep is activated, so at low temperatures/stresses; it equals 1 when only dislocation-creep acts, so at higher temperatures/stresses.
The total dissipative potential for stationary creep j vp ðs; TÞ is then expressed as follows: where subscript 1 stands for the linear cation diffusion creep and 2 for the non linear dislocation creep. The description of mechanical behavior of the equivalent solid material is completed by accounting for the elastic strain. Assuming after Leclerc [17] the additivity of macroscopic strains under the hypothesis of small strains for nuclear pellet, the macroscopic mechanical stresss for the solid homogenized material can be expressed as: where C hom is the effective elastic tensor of the equivalent homogeneous material taking into account the presence of porosity,ε is the macroscopic strain tensor, which is the result of an elastic contributionε e and an inelastic one, given by the viscoplastic straiñ ε vp . Martin [18] initially derived correlations to express material properties as a function of the temperature, which were then used in 3D numerical simulations carried out by Gatt [19] to derive the equivalent Young's modulus and shear modulus as functions of both temperature and porosity fraction Eðf ; TÞ and Gðf ; TÞ for T 2929 and porosity f 0.3. The expressions reported as follows shall be used in section 5: where ða i ; b i ; c i Þ are material parameters, experimentally defined.
3. Thermal-physical characterization of the melting process

Parameters influencing the melting interval
Fuel material is an heterogeneous composition of a solid and a liquid phase during melting. The solid-liquid phase transition occurs over a melting interval characterized by a lower melting temperature, the ''solidus temperature" T sol , and an upper melting temperature, the ''liquidus temperature" T liq .
According to literature [20], the melting interval of nuclear fuel material essentially depends on three parameters: i. In the specific case of mixed oxides of uranium and plutonium (the so called MOX fuel) the composition of the material at fabrication in terms of volume fraction of PuO 2 and UO 2 ; ii. The burnup, which describes the amount of uranium burnt during irradiation and it is expressed in GWd/tU; iii. The ratio oxygen/metallic actinides (O/M) present in the fuel.
We detail in the following how each of those parameters impact the melting process of non stoichiometric uranium dioxide imposing a modification of the solidus or liquidus temperatures which are given by: where T f is the melting temperature of stoichiometric uranium dioxide. This temperature has been recently re-evaluated by Manara [21] who confirmed previous measurements and proposed an updated melting temperature for the stoichiometric UO 2 which amounts to 3147 K ± 20 K, valid in the pressure range from 10 to 250 MPa.

The plutonium content
Depending on the fuel material type, the volume fraction of plutonium in the nuclear fuel can evolve considerably during irradiation, from 0% up to 25%. The correlation describing the evolution of melting interval of (U,Pu)O 2 material with the composition dates back to experimental measurements carried out in the 70's by Leon and Bailly [22] and by Aitken and Evans [23]. Melting temperatures of stoichiometric plutonium dioxide have been recently re-evaluated by De Bruycker [24]. He highlighted a difference with previous experimental characterizations of approximately 300 K. De Bruyker determined a melting temperature amounting to 3017 K ± 20 K while previous characterizations gave a melting temperature of 2701 K±35 K, according to Carbajo and al. [25], and 2660 K, Gu eneau and al [26]. In its article, De Bruycker details that the difference is due to the previous experimental techniques, which are affected by extensive interaction between the PuO 2 samples and the tungsten containment.The recommended solidus and liquidus curves for a stoichiometric fresh fuel containing a certain volume fraction of Plutonium, y, result from detailed analysis [27] of experimental data [22,23] obtained previous to De Bruycker study and reported in Ref. [ (9) where y corresponds to the volume fraction of plutonium.

Burn-up
The effect of burnup on the melting interval has been largely investigated in past years by Ref.
[27e30] and recently by Carbajo [25]. According to Carbajo review of Bu effect on the melting point of stoichiometric UO 2 and in agreement with the evaluation proposed by Oak Ridge National Laboratory [20], a decrease of T sol of 0.5 K each GWd/tU is recommended:

Effect of oxygen to actinides ratio
According to past studies [21,30] and more recent analysis carried by Gu eneau [31], the evolution of T sol and T liq as a function of O/ M are showed in Fig. 2. The oxygen to actinides ratio O/M of the fuel material is a highly influencing parameter: the deviation from stoichiometry induces a decrease in T sol and T liq . Also it progressively increases the difference between the two, thus leading the melting interval to occur over a wider range of temperatures.
During irradiation, for each fission of one atom of uranium, two atoms of oxygen are freed and combined with some Fission Products (FPs). The majority of FPs are anyway noble gasses and rare ones which are characterized by a low probability of combining with oxygen. As a consequence the O/M ratio increases during the life of the fuel material under neutronic flux. Studies carried out by Riglet-Martial et al. [32] showed that thermal gradient characterizing the fuel material during irradiation is the source of a strong Evolution of T sol and T liq as a function of the oxygen to actinides ratio in the UO 2 nuclear fuel. Data from Ref. [21] for hyperstoichiometric UO 2 and from Ref. [30] for hypostoichiometric. The deviation from stoichiometry reduces the melting temperature and progressively increases the melting interval.
reduction-oxidation reaction leading to migration of the oxygen from the center of fuel pellet to its periphery which induces a decrease of O/M ratio in the central zone. It has been showed by an expert group and then reported by the French Autority of Nuclear Safety [33] that in a hypothetical accidental case leading to fuel melting, this shall appear in the central region of the PWR fuel pellet and progressively evolve towards the periphery. Melting shall thus take place in a hypostoichiometric system.
For the hypostoichiometric region of UO 2 , the recommended correlations in the validity region 1:94 O=M 2:00 are the following [

The melting interval
According to previous discussion, we expect the central melting region to have an hypostoichiometric composition (O/M < 2.00). The correlations to describe the evolution of solidus and liquidus temperatures as a function of the composition and of material irradiation have been proposed by Refs. [30]. Under the hypothesis that the three effects can be linearly superposed, the melting interval can be expressed as follows: where all the symbols have already been defined.

Discontinuity in material density
Solid-liquid phase transition of UO 2 is characterized by a discontinuity in material density. Recently, Fink [34] reviewed all correlations concerned with UO 2 . For the evaluation of liquid UO 2 density, the recommended correlation is the one proposed by Breitung and Reil [35] which turns out a density of liquid UO 2 at melting temperature of 8860 kg:m À3 . For the density of solid UO 2 he recommended the correlation obtained by Martin [18], valid for temperature 923 K T 3120 K. According to this evaluation, density at melting temperature amounts to 9560 kg:m À3 ± 40 kg:m À3 and the discontinuity at melting point for a stoichiometric uranium dioxide amounts to 700 kg:m À3 ± 40 kg:m À3 . The discontinuity in material density turns into a volume expansion DV V m obtained as follows: where Dr ¼ r sol À r liq is the difference in density between the solid and liquid UO 2 at melting temperature and r m ¼ r sol . The characteristic expansion occurring at melting and associated to the discontinuity in material density, amounts to about 7.32% ± 0.39% according to Fink recommendation.

Methodology and hypothesis
We consider the onset of melting in a PWR type fuel pellet. Due to the thermal profile that establishes during irradiation, the highest temperatures are reached always in the pellet center. At a microscopic scale, the liquid inclusions originate in the pellet center when the temperature locally overcomes the melting threshold of T sol . Under the assumption that they are randomly distributed and of spherical shape, the melting material can be assumed to be isotropic and the liquid presence can be completely described in terms of its volume fraction, x liq .
The thermal-mechanical behavior of the solid has been exhaustively described in the previous sections: we consider its elastic properties as described in equation (6) to remain constant and equal to the value obtained at T sol during the melting process. One should note that the validity range of those correlations is limited to 2929 K. No correlations are available for the temperature interval [2929e3147] K, being 3147 K the measured melting temperature of stoichiometric UO 2 . We assume that no unreasonable values of elastic properties should be obtained in this range by means of previously proposed relations.
In our study we consider the forming liquid material as incompressible newtonian liquid: the bulk modulus of the liquid UO 2 K l /∞. We can easily verify that the hypothesis that the shear modulus of liquid UO 2 G l /0 is coeherent with material properties at melting points. If we evaluate the dynamic viscosity of uranium dioxide at melting temperature by the equation proposed by Woodley and presented by Fink in Refs. [34], we find it amounts to 4.32 cP (mPa.s) and thus negligible with respect to solid one at solidus temperature, which is in the order of the GPa. We report for completeness the mentioned correlation for the dynamic viscosity m l [34]: where T f is the melting temperature of stoichiometric UO 2 , 3147 K. Liquid fraction xliq is expected to vary according to the composition of material undergoing melting, which is temperature dependent. For the generic temperature T sol T T liq during phase transition interval, the generic composition c of the material can be expressed in terms of the solidus and liquidus compositions c sol and c liq respectively, weighted on the volume fraction of each phase: We propose to derive the liquid fraction expression from correlations of liquidus and solidus temperatures through the lever rule: for the generic temperature T in the melting interval the liquid volume fraction x liq is function of solidus and liquidus temperatures previously expressed in equations (12) and (13).
In our study we consider the PWR type fuel, where the volume fraction of plutonium is limited to a negligible amount so that we can assume that DT y sol ¼ 0 and DT y liq ¼ 0. We apply the lever rule in the thermal interval T sol T T liq where the liquid fraction can be expressed as a function of generic stoichiometry O/M and burnup Bu x liq ðT; BU; O =MÞ: x liq ðT; BU; O = MÞ ¼ 10 3 ðO=M À 2:00Þ À 0:5BU þ T f À T

2:57
T À T f À 0:5Bu Fig. 3 illustrates the effect of the burnup and of the stoichiometry on the melting interval. In (a) the evolution of the liquid fraction with the temperature in the melting interval is evaluated for a fixed value of stoichiometry and for three values of burnup: 0 GWd/tU corresponding to a non irradiated fuel, 30 GWd/tU and 50 GWd/tU. We observe that the melting process occurs over a thermal interval because the T sol is reduced by the rise of the burnup. In figure (b) of 3 is illustrated the effect of different stoichiometries on the melting interval, for a fixed burnup. The decrease of stoichiometry results into a reduction of the solidus threshold while the thermal interval is not affected.

The RVE, the Representative Volume Element
We apply the homogenization approach to analytically evaluate macroscopic behavior of heterogeneous materials. It consists in replacing the real heterogeneous one with an homogenized medium which behavior is macroscopically equivalent to the real one. We define a Representative Volume Element (RVE), denoted V , which is large enough to consider the melting UO 2 nuclear fuel as an isotropic matrix including the liquid portion, being characterized just by its volume fraction as defined in eq (17). According to the theory of effective moduli as described by Gilormini and Bornert in Ref. [36] we can define the macroscopic stress fields and the macroscopic strain fieldε obtained by averaging on the RVE V the microscopic stress tensor s and the microscopic strain tensor ε which are space dependent: Wheres andε account for the presence of liquid in the material. Under the assumption of small strains and according to the approach used by Monerie in Ref. [6], we assume the additivity of macroscopic strains to stand. The thermal-mechanical strain for the equivalent melting material can be thus expressed as the sum of the elastic strainε e , the viscoplastic strainε vp , the thermal dilatationε th . Also, we need to account for the characteristic strain due to density reduction, defined in the followingε m . The total strainε is:ε ¼ε e þε vp þε th þε m The equivalent homogenized stress tensor can be thus expressed as: where C hom ¼ ð3K; 2GÞ is the effective elastic tensor,K andG are the effective bulk modulus and the effective shear modulus of the equivalent material, respectively.

The thermal-elastic behavior
According to classic writing the thermodynamic potential J of the equivalent material can be written as follows: Where r is the material density,ε H ¼ 1 3 trðεÞ is the hydrostatic component of the strain tensor,ẽ ¼ε Àε H I is the deviatoric component of the strain tensor and I the 4 th order identity tensor,ã is the macroscopic dilatation coefficient for the equivalent homogenized material,c p is the equivalent heat capacity of the homogenized material, T 0 is the reference temperature.In addition to the usual terms of elastic deformation and thermal dilatation, during fuel melting due to the discontinuity in material density, the thermal-elastic strain presents an additional contribution due to the characteristic decrease in density at meltingε m . Under the assumption of isotropic strain, we can express this term under the form of an equivalent thermal dilatation coefficient which will be further detailed in section 4.3.3: The thermodynamic potential J is modified as follows: The macroscopic mechanical behavior of the equivalent material s can be obtained by deriving the potential with respect to the strain tensor: The specific entropy of the system results from the derivative of the potential with respect to the temperature T. Taking into account the dependency of material properties with respect to temperature T we have: In other studies, where it is generally considered that the derivative with respect to time of material properties are negligible, this assumption leads to obtain the usual form of the specific entropy: wherec p is the equivalent heat capacity of the melting material, defined as a function of specific heat capacity of the solid at solidus temperature c p;s ðT sol Þ, liquid fraction and latent heat L: In the generic case we can indeed consider the thermal balance for the unitary melting mass in a closed system: where hðtÞ is the enthalpy of melting material at time t during transition process and q the specific heat actually inducing the melting process. By taking advantage of the extensive nature of enthalpy we can say that: where h sol is the enthalpy of the material at solidus temperature and L the latent heat as previously mentioned. Thus replacing in eq. (30) and taking advantage of the relation between enthalpy and temperature through specific heat capacity we can easily obtain eq. (29).

The equivalent elastic strain
According to previously stated hypothesis, the effective elastic tensor of the melting fuel material C hom ¼ ð3K; 2GÞ is obtained considering a newtonian incompressible liquid phase having K l / ∞ and G l /0 and a solid phase whose properties are obtained according to eq. (6), evaluated at T sol . Fig. 4 shows the evolution of effective properties for dense UO 2 . Graphs are normalized with respect to the solid property for the shear modulus and with respect to the liquid property for the bulk modulus. They are obtained according to different approaches: Voigt and Reuss upper and lower bound respectively are evaluated to identify the region in which we expect the evaluation to vary, then the Hashin and Shtrikman (HS) approach is applied according to the Mori-Tanaka (MT) and Self-Consistent (SC) schemes [36]. The general formula proposed by HS for the evaluation of bulk and shear modulus is proposed as follows: where p stands for the generic phase, solid or liquid, C:D indicates the average, K * and G * refer to the reference medium and are defined as . To apply the MT scheme, the reference medium corresponds to the solid matrix, whereas the SC is based on the hypothesis that the reference medium is the homogenized system itself.
As showed in figure (a) of 4, the MT and SC estimations are coherent as long as liquid fraction xliq < 0:1. The self-constistent scheme based on the estimation of Hashin and Shtrikman as reported in Ref. [36] is retained for the evaluation of the effective properties as it is better adapted for higher liquid fractions. Once solved eqs. (34) and (35) the elastic properties for the equivalent elastic material are the following: We observe that vG=vx liq 0 and is null for x liq ¼ 0:6, which is the well known solution for an heterogeneous material having incompressible matrix and an internal incompressible void fraction [36]. Even though the bulk modulus tildeK following homogenization depends on the equivalent shear modulus, the fact that K s [ G s induces no relevant dependencies. If we had a different case, we would have a non monotonic function of the liquid fraction as we immediately see once considering its derivative with respect to x liq : Thus we have a maximum of the bulk modulus for x liq < 0:6 and for x liq > 0:6 it increases as the square of x liq .

The equivalent thermal strain
Under the hypothesis of isotropic mechanism, the thermal dilatation impacts the hydrostatic behavior of the medium: According to Levin relation as reported in Ref. [37], the equivalent dilatation coefficient of an heterogeneous material can be expressed as a function of the equivalent bulk modulusK:ã where CaD ¼ ð1 Àx liq Þa s þx liq a l is the average dilatation coefficient.

The characteristic melting expansion
Still assuming the occurrence of an isotropic phenomenon, the characteristic strain due to density reduction impacts the hydrostatic term of the macroscopic behavior ε m ¼ε m I Volume expansion associated to discontinuity in density concerns just the melted volume fraction. The linear strain can be expressed as:

The viscoplastic mechanical behavior for an incompressible matrix
As presented in section 2, the total viscoplastic behavior of the system is given by the diffusion and the dislocation mechanism, modeled and coupled according to eqs. (4) and (3). The non linear strain rate for viscous solid material is obtained by deriving eq. (4) with respect to the macroscopic stress tensor [6]: We assume as a first order approximation that, due to the high temperatures reached by the fuel at melting onset, porosity in the central region can be neglected. This results in imposing A i ðf Þ ¼ 0 in eq. (3). Under this assumption the i th creep mechanism can be expressed according to a Norton potential. The resulting strain rate of the solid matrix at melting onset is in the form: Hðn i Þs n i eq s s eq (45) Before evaluating the behavior of the equivalent homogeneous melting material, linearization of the creep law is obtained by applying the modified secant method as presented by P. Suquet in Ref. [36]. The strain can be locally described for each point x2 V as a function of the local stress: εðxÞ ¼ L sct ðx; sðxÞÞ : sðxÞ (46) where L sct ðx; sðxÞÞ ¼ C sctÀ1 ðx; εðxÞÞ is the effective secant tensor which depends on the local strain εðxÞ and on the local phase (solid or liquid). We evaluate the global non linear behavior on the basis of the approximate local solution: average strain and stress in the non linear heterogeneous material are considered equal to those of the linear equivalent medium. The macroscopic global strain is thus: whereLðc r ; L sct r ; …Þ is the effective tensor of the linear equivalent medium. It is obtained through the secant model on the constitutive model of the non linear material and on the volume fraction of components.
According to the modified extension of the secant method, we identify an effective secant tensor for each phase r as a function of the second order moment of stress tensor.
J and K are the 4 th order tensor defined as K ¼ IÀ J where I is the 4 th order identity tensor and J abgd ¼ 1 2 i ab i gd . When discussing the equivalent viscoplastic behavior of the homogenized medium, we evaluate the particular case in which we have a single non linear phase (the solid matrix) and liquid inclusions which can be assumed to be isotropically distributed in the solid. It first results that the only variable to be evaluated is the secant shear modulus of the non linear solid phase. This is defined through the non linear behavior of the solid phase. Thus after eq. (45) we have: Then, the linear equivalent material can be characterized by effective properties (3K; 2G) obtained considering that both components are incompressible and that shear modulus of the liquid G l /0: Then, the second order moment of the equivalent strain and stress tensors were obtained by P. Suquet in Ref. [36], simplified as both phases are incompressible: wheres eq is the macroscopic equivalent stress. Solving eq. (51) with respect to the effective shear modulus and secant one defined in eq. (50) and eq. (49) respectively, we obtain the model describing the mechanical behavior of the equivalent medium: describing the effect on the mechanical behavior due to the presence of the liquid fraction x liq . Parameters for the two mechanisms have been discussed by Monerie and Gatt in Ref. [6]: we retain that n 1 ¼ 1 for scattering mechanism and n 2 ¼ 8 for dislocation mechanism. Fig. 6 shows that bðx liq ; n i Þ acts like a multiplier of strain rate and that dislocation-creep is the driving phenomenon leading to a rapid loss of mechanical resistance of the melting volume.

Definition of the test case
We test the model in a simple case in order to easily access and study the calculated behavior. For this we use a 2D axis-symmetric system having the geometrical characteristics of the fuel pellet. The material properties are the ones of irradiated UO 2 , with stoichiometry O/M ¼ 1.98 and Bu ¼ 30 GWd/tU and null plutonium content. The problem is solved by means of the finite element code Cast3m while material properties are treated via MFRONT. Fig. 7 presents geometry and boundary conditions of the problem: L0 corresponds to the central axis; LB and LH are respectively the bottom and top edge, free to move only in the radial direction. Their length corresponds to pellet radius. LR is the external edge of the pellet, which length equals pellet height.
The system is submitted to a parabolic radial temperature profile which ensures that melting first origins in the central zone and then progressively moves towards pellet periphery. At the periphery of the system the temperature is fixed. We imposed a linear rise of the thermal profile according to a constant rate.

Convergence test: node density
The system is modeled with a uniform radial density nDr with a ratio Dr=Rz3:10 À3 . The reduced radial size is required as a consequence of the strong thermal gradient, to have a minimum of 2 nodes that undergo melting. Fig. 8 shows the convergence of simulation outputs as a function of the radial mesh density. Errors are calculated for the edge displacement and the equivalent stress tensor at r ¼ R and T ¼ 3200 K with respect to the highest tested density (n ¼ 6). The model converges with a residual error in the order of 10 À4 for the displacement and 10 À3 for stresses.

Evolution of radial stress
We have presented in the article the evaluation of the elastic and viscoplastic behavior. Fig. 9 illustrates the different evolution of the equivalent stress tensor s eq along the radial dimension. One the left it is evaluated according to a pure elastic behavior and to an elasto-viscoplastic one, on the right. Both are evaluated at 3200 K. The highlighted region corresponds to the melting plateau that, according to the defined thermophysical characteristics of the material, takes place between 3112 and 3141 K.
As the liquid material is modeled as incompressible, s eq is null in the melted region. In the solid region, we observe that the creep mechanisms activated at high temperature induce a strong relaxation of stresses. Toward the external side, temperature decrease and stresses increase in the system due to the absence of creep.

The resulting external displacement
We are interested in evaluating how the macroscopic displacement of the system varies according to the central melted region. Fig. 10 shows the kinetics of the modeled system edge displacement as a function of time. Dt is defined as the difference between the generic time t within the transient and the t 0 at which liquidus temperature is first reached in pellet center.
The melted fraction is defined as the ratio of melted volume with respect to the total volume of the pellet: where r liq is the radius at which the liquidus temperature is reached. Calculation shows that once melting appears, a discontinuity in displacement kinetics is observed. This leads to an acceleration of system deformation. In the proposed problem we imposed a linear increase of temperature equal to 1 K=s and the kinetics of displacement is accelerated by a factor 3.
The implemented model shows that there is a direct relation between the calculated displacement kinetics and the increase of the melted fraction, which is the major driver to its evolution. From the figure we can conclude that the thermal-mechanical strain defined in equation (21) can be approximated by _ εz _ ε m .

Conclusions and further work
In this article we develop and numerically implement a model to describe the thermal-mechanical behavior of the UO 2 fuel material undergoing initial central melting. The purpose is to validate the modeling principle and study the calculated modifications in the macroscopic behavior. We start from the viscous mechanical behavior of the solid porous fuel material obtained by elliptical thermal-mechanical potential as described by Monnerie and Gatt in Ref. [6]. Fundamental to our model is the assumption that porosity is negligible as the fuel approaches the melting temperature. Thus, the viscoplastic strain rate for the solid fuel material at melting onset is obtained from a Norton potential. Liquid inclusions progressively nucleate in the central region and their effect is evaluated by means of homogenization approach. Under the assumption that liquid inclusions in the volume undergoing melting appear according to a uniform distribution, the liquid phase is completely described by the volume fraction x liq . We Fig. 7. Definition of the thermal-mechanical test case: geometry (left) and boundary conditions and thermal-loading (right). Fig. 8. Convergence of the displacement u(r) and of the equivalent stress tensor evaluated at r ¼ R and t ¼ 3200 K, with respect to the mesh density. evaluate the macroscopic elastic properties of the equivalent material as a function of liquid fraction, using the Self-Consistent model. We account for material density assuming that an isotropic expansion originates and thus introducingε m contribution to the macroscopic strain rate. Concerning the viscoplastic behavior, we use modified secant method to obtain the macroscopic behavior and we demonstrate that the viscoplastic deformation undergoes strong acceleration due to dislocation creep, leading to a strong decrease in system stress. The model is implemented in Cast3m and the macroscopic mechanical behavior is analyzed in a simple test case. We study a 2D axisymmetrical system having the radius and height of the fuel pellet, blocked at the top and bottom and free to deform in the radial direction. The modeled material properties are the ones discussed for the UO 2 material and are temperature dependent. By comparing a pure elastic solution and the elasto-viscoplastic one we show that the creep mechanisms create a strong relaxation of stresses with respect to the elastic case, with an increase of calculated stresses only in the external part of the pellet, as a consequence of the decrease of temperatures. Through the analysis of the macroscopic displacement it is evident that the discontinuity appearing at melting onset and the resulting kinetics are directly linked to the appearance and progression of melting in the pellet center.
The current work is dedicated to the implementation of the model in Alcyone thermal-mechanical code, where the characteristic evolution of stresses associated to the melting process have to be coupled with the characteristic phenomena occurring in the fuel pellet, like swelling and relocation. Further work is required but it takes advantage of the fact that the melting model is developed according to the same approaches and using the same independent variables. The modeled mechanisms of the fuel element will result from the tendency of the melted fraction to expand, as we have shown in Fig. 10 but also on the power history of the rod, controlling the available latent heat to melt, on the evolution of gaseous swelling at high temperature and on the reaction of the surrounding cladding to the imposed stresses.

Declaration of competing interests
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.   Parrat: Writing -review & editing, Supervision.