Depletion, Chemical Reaction and Transport in High Burnup Nuclear Fuel

ORNL/SPR-2019/1104 Depletion, Chemical Reaction and Transport in High Burnup Nuclear Fuel S. Simunovic J. W. McMurray T. M. Besmann E. E. Moore K. T. Clarno W. A. Wieselquist M. H. A. Piro Approved for public release. Distribution is unlimited. March 2019 DOCUMENT AVAILABILITY Reports produced after January 1, 1996, are generally available free via US Department of Energy (DOE) SciTech Connect. Website http://www.osti.gov/scitech/ Reports produced before January 1, 1996, may be purchased by members of the public from the following source: National Technical Information Service 5285 Port Royal Road Springfield, VA 22161 Telephone 703-605-6000 (1-800-553-6847) TDD 703-487-4639 Fax 703-605-6900 E-mail [email protected] Website http://www.ntis.gov/help/ordermethods.aspx Reports are available to DOE employees, DOE contractors, Energy Technology Data Exchange representatives, and International Nuclear Information System representatives from the following source: Office of Scientific and Technical Information PO Box 62 Oak Ridge, TN 37831 Telephone 865-576-8401 Fax 865-576-5728 E-mail [email protected] Website http://www.osti.gov/contact.html This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. ORNL/SPR-2019/1104 Computational Sciences and Engineering Division DEPLETION, CHEMICAL REACTION AND TRANSPORT IN HIGH BURNUP NUCLEAR FUEL S. Simunovic1 T. M. Besmann2 E. Moore2* K. T. Clarno1 W. A. Wieselquist1 J. W. McMurray1 M. H. A. Piro3 1Oak Ridge National Laboratory, Oak Ridge, TN of South Carolina, Columbia, SC 3University of Ontario Institute of Technology, Oshawa, ON, Canada *Currently at Lawrence Livermore National Laboratory 2University Date Published: March 31, 2019 Prepared by OAK RIDGE NATIONAL LABORATORY Oak Ridge, TN 37831-6283 managed by UT-BATTELLE, LLC for the US DEPARTMENT OF ENERGY under contract DE-AC05-00OR22725 CONTENTS CONTENTS.................................................................................................................................................III ABSTRACT...................................................................................................................................................1 1. INTRODUCTION .................................................................................................................................1 2. MASS AND HEAT TRANSPORT IN CHEMICAL AND TEMPERATURE GRADIENTS.............3 2.1 TRANSPORT MODEL BASED ON IRREVERSIBLE THERMODYNAMICS .....................3 2.2 MECHANISTIC MODEL FOR MASS DIFFUSION ................................................................5 2.2.1 Mass flux representation ................................................................................................5 2.2.2 Conservation equation for mass transport......................................................................6 3. FEM IMPLEMENTATION OF THE MECHANICSTIC TRANSPORT EQUATIONS.....................7 4. MODELING OF THE FUEL COMPOSITION WITH BURNUP .....................................................10 5. IMPLEMENTATION OF TRANSPORT MODEL IN BISON ..........................................................10 6. EXAMPLES.........................................................................................................................................11 6.1 OXYGEN TRANSPORT DRIVEN BY CHEMICAL POTENTIAL GRADIENT .................12 6.2 OXYGEN TRANSPORT DRIVEN BY CONCENTRATION GRADIENT...........................16 7. SUMMARY .........................................................................................................................................18 8. REFERENCES.....................................................................................................................................18 iii ABSTRACT We have developed a formulation and computational model for oxygen transport in Light Water Reactor (LWR) uranium dioxide fuel. The overall model couples the burnup simulation isotopic composition with a thermochemistry model of the fuel phase, and oxygen transport using the driving forces from the thermodynamic calculations. The diffused oxygen is accounted for in the thermochemistry model in order to establish consistent material and thermodynamic conditions in the fuel undergoing burnup. The model has been implemented in the nuclear fuel performance code Bison. The formulation and the model have been successfully demonstrated on fuel burnup data from the open literature. The developed capability enables consideration of complex chemical composition of irradiated fuels beyond available burnup models in Bison. 1. INTRODUCTION The oxidation state and evolving composition of light water reactor (LWR) nuclear fuels have a strong influence on the performance and safety of fuel elements and thus reactor operation (Olander, 1976). They affect almost all of the processes of practical importance. The temperature and the oxygen chemical potential determine the chemical form of the fission products, e.g., whether they are stable as a metal or an oxide, dissolve in the fuel matrix or form new phases, etc. Characterization of the fuel composition during burnup have been presented in numerous publications as summarized in (Olander, 1998, Walker, Rondinella, Papaioannou, Van Winckel, Goll and Manzel, 2005, Park, Yang and Park, 1997) and will not be reviewed in this document. The importance of the modeling of fuel oxygen transport under large composition and temperature gradients (Marchant and Bowen, 1975, Adamson, Aitken, Evans and Davies, 1975, Adamson and Carney, 1974, Evans, Aitken and Craig, 1969, Marin and Contamin, 1969, Janek and Timm, 1998, Lassmann, 1987, Olander, 1972, Matzke, 1987) is taken here for granted and the main focus is on describing the model formulation and implementation. Models for fuel composition are usually derived from experiments on fresh or simulated irradiated fuel, which may rely on the assumption of a single UO2 solid solution phase (Lindemer and Besmann, 1985). In reality, as the irradiation process progresses, multiple phases may form due to the creation of fission products from the consumed actinides and their interactions with the fuel matrix phase. The models based on the thermodynamic equilibrium yield the phase assemblage and phase compositions that represent the lowest free energy for a given temperature and pressure (Besmann, McMurray and Simunovic, 2016, Gueneau, Baichi, Labroche, Chatillon and Sundman, 2002, Gueneau, Dupin, Sundman, Martial, Dumas, Gosse, Chatain, De Bruycker, Manara and Konings, 2011). These models are assumed to be a good approximation of the real phenomena because the high temperature of the fuel and long-time intervals at steady power generation are driving the fuel to approach thermodynamic equilibrium, at least locally (Besmann, 2012). Thermodynamic equilibrium models by definition do not account for the kinetics of chemical processes that lead to the equilibrium state. However, they can provide the driving forces for kinetic models (Demirel, 2007, Tschoegl, 2000), such as those for oxygen transport. We have developed a thermodynamic computation library, Thermochimica (Piro, Simunovic, Besmann, Lewis and Thompson, 2013), for modeling composition and thermodynamic properties of complex material systems. It accommodates materials with large numbers of components and phases that must 1 typically be considered in irradiated nuclear fuel. Figure 1 shows an example of a computation of the radial distribution of phases across a nuclear fuel pellet (Piro, Banfield, Clarno, Simunovic, Besmann, Lewis and Thompson, 2013) at a burnup of 102 GW·d·t(U)-1 reflecting the experimental observations from (Walker, Rondinella, Papaioannou, Van Winckel, Goll and Manzel, 2005). Figure 1. Predicted radial distribution of phases across an LWR fuel pellet at an average burnup of 102 GW·d·t(U)-1. As is to be expected, the fluorite oxide phase is dominant with several additional minor phases predicted to be stable, including the observed noble metal inclusions in the face centered cubic and hexagonal closed packed crystal structures, plus secondary oxides, and vapor species. The co-existence of these phases is important in computing the oxygen-to-metal ratio and oxygen chemical potential. The fission, radioactive decay products, and transuranics together with high temperatures result in the formation of structural defects in the primary fuel phase. These defects in the UO2 lattice control a wide range of phenomena in the nuclear fuel. The dominating structural disorder created by irradiation and temperature is the formation of Frenkel defects on the anion, oxygen, sub-lattice (Matzke, 1981). Schottky defects are less prevalent and occur on both anion and cation sub-lattices. Oxygen transport occurs primarily through coordinated movement of defects on the anion lattice, whereas the movement on the cation (i.e., uranium) lattice is much slower but significant enough to control the creep rate by movement of uranium vacancies. Of primary interest here is the development of models for oxygen transport within the fluorite phase in large compositional and thermal gradients. In such setting, the diffusion is governed by more than just a gradient of the concentration. The diffusing species chemically interacts with the transport matrix and is also influenced by variations in temperature and pressure. For the compositions and temperatures of interest the fuel is essentially single-phase fluorite structure dioxide, although minor noble metal phases are observed at high temperature, no other oxide phases are observed with the possible exception at the very periphery of a fuel pellet. Thus, oxygen diffusion through only the single-phase material need to be 2 considered. In an earlier report, we have derived the oxygen transport model based on irreversible thermodynamics (de Groot and Mazur, 2011). Here, we take a more mechanistical approach that is based on the minimization of Gibbs energy and the driving forces from the thermodynamic equilibrium model (Kocherginsky and Gruebele, 2016). The recent developments in modeling diffusion and continuum representations of irreversible thermodynamics in control volumes were used to formulate the computational model for transport in concurrently varying composition and temperature. Models for the mobility of oxygen in hypo- and hyper-stoichiometric nuclear fuel from recent work by Moore et al. (Moore, Gueneau and Crocombette, 2013) were used. The changing fuel material elemental composition during burnup was calculated using the Origen software. The simulated composition was read into the oxygen transport model using the solution function method in the fuel performance code Bison (Hales, Novascone, Spencer, Williamson, Pastore and Perez, 2014). Our transport model implementation in Bison is computationally expensive because the thermodynamic equilibrium calculations involve a nested iterative optimization at every computation point in space. Its practical deployment will require studying the rate of compositional changes and their effects on other fuel models. We are also developing new algorithms for improving the solver performance in the Finite Element Method (FEM) (Zienkiewicz, Taylor and Zhu, 2013) simulation framework which will make Thermochimica more computationally affordable for other aspects of fuel performance modeling. 2. 2.1 MASS AND HEAT TRANSPORT IN CHEMICAL AND TEMPERATURE GRADIENTS TRANSPORT MODEL BASED ON IRREVERSIBLE THERMODYNAMICS Thermodynamic models calculate the composition, chemical, and physical properties of a material system that is in equilibrium (Hillert, 2007). The Gibbs energy is at a minimum at constant temperature and pressure, with no entropy production and no thermodynamic forces operating within the system. However, heat and mass transport are non-equilibrium dissipative processes, driven by the unbalanced driving forces. The processes considered here are assumed to occur at sufficiently large time intervals to be treated using steady-state irreversible thermodynamics (Tschoegl, 2000). The systems are open and thus can exchange energy and mass with the environment. Constant driving forces result in stationary fluxes and stationary states. For example, the temperature or composition profile does not change under steady-state conditions, although entropy is being produced in the system. The constant driving forces and fluxes result in a steady state irreversible system, which can be described with equilibrium thermodynamics models. In effect, equilibrium thermodynamics is used to model kinetic, non-equilibrium phenomena (Kondepudi and Prigogine, 2015), provided steady-state is achieved. It is implicitly assumed that the system is sufficiently near equilibrium that there is a linear relation between the driving forces and rate processes. These assumptions are less restrictive in transport models which assume that the gradients are not so large that they influence the chemical reaction models. Heat and mass transport simulations are usually cast in the framework of irreversible thermodynamics. The method is based on the assumption of local equilibrium in the constitutive volumes, and conservation of energy and mass in the system. The macro non-equilibrium system is assumed to be an assembly of open elemental volumes each at equilibrium, and thus equilibrium thermodynamic relations are valid for locally defined thermodynamic variables. The elemental volumes can be used to describe heterogeneous systems as long as the temperature can be well-defined at every location. The intensive thermodynamic 3 variables temperature (T), pressure, (P) and chemical potential of substance k (𝜇𝑘), become functions of position x, and time, t: 𝑇 = 𝑇(𝒙,𝑡) , 𝑃 = 𝑃(𝒙,𝑡) , 𝑢 = 𝑢(𝒙,𝑡) , 𝜇𝑘 = 𝜇𝑘(𝒙,𝑡) (1) The extensive thermodynamic variables are replaced by their volumetric densities as: 𝑠 = 𝑠(𝒙,𝑡) , 𝑛𝑘 = 𝑛𝑘(𝒙,𝑡) (2) where s, u, and 𝑛𝑘 denote entropy per unit volume, internal energy per unit volume and moles per unit volume of substance k, respectively. Integrated values of the densities over the system volume do not mutually correlate by standard thermodynamic equilibrium relations because the system is not in equilibrium. However, the local thermodynamic equilibrium relations are valid as long as temperature and composition are well defined for each point in space and time. The nonequilibrium thermodynamic formulation (Tschoegl, 2000) is cast in the form of internal entropy production density, 𝜎: 1 ― 𝑇 () 𝜎 = 𝐽𝑈 ∙ ∇ ∑ 𝑘 𝜇𝑘 () 𝐽𝑘 ∙ ∇ 𝑇 (3) where 𝐽𝑈, and 𝐽𝑘 denote flux of energy and flux of substance k, respectively. The energy flux and temperature distribution are known, and the mass flux is determined by the continuity equation: ∂𝑛𝑘 ∂𝑥 + ∇ ∙ 𝐽𝑘 = 0 (4) In linear irreversible thermodynamics, the fluxes are linearly proportional to the driving forces and the proportionality factors are phenomenological expressions that do not depend on the gradient values of the thermodynamic variables. In the linear regime, the system evolves to stationary, steady state with a constant entropy production. When fluxes are generalized to include coupling terms of the same dimensionality, proportionality coefficients 𝐿𝑖𝑗 couple all the driving forces: 𝐽𝑖 = ∑𝐿 𝐹 = 𝐋 𝐅 (5) 𝑖𝑗 𝑗 which form an Onsager (Onsager, 1931) matrix, L, and a vector of driving forces, F. Assuming only one specie k is transported, its mass flux is then defined: 𝜇𝑘 () 𝐽𝑘 = ― 𝐿𝑘𝑘∇ 𝑇 1 𝑇 () ―𝐿𝑘𝑞∇ The proportionality coefficients for mass transport equation can be written as products of the species concentration and the kinetic coefficient: 4 (6) 𝜇𝑘 () 𝐽𝑘 = ― 𝑛𝑘𝐿𝑘𝑘∇ 1 𝑇 () ― 𝑛𝑘𝐿𝑘𝑞∇ 𝑇 (7) In order to separate the temperature from the chemical potential gradient term in Eq. 7: 𝜇𝑘 𝜇𝑘 1 = ∇𝜇𝑘 ― 2 ∇𝑇 𝑇 𝑇 𝑇 () ∇ (8) where the gradients on the right-hand side are not restricted to constant temperature or composition. Gradient of the chemical potential is usually expressed as a sum of partial gradients with respect to composition and temperature which leads to: 𝐽𝑘 = ― 𝑛𝑘𝐿𝑘𝑘 ∂𝜇𝑘 𝐿𝑘𝑞 ∂𝜇𝑘 1 ∇𝑛𝑘 + ― 𝜇𝑘 ― 𝑇 ∇𝑇 𝑇 ∂𝑥 ∂𝑇 𝑇𝐿𝑘𝑘 𝑇 { [ ( )] } (9) The first term in curly braces of Eq. (9) represents the customary diffusion term due to a concentration gradient. The expression in the square braces describes the effect of a thermal gradient. ∂𝜇𝑘 While the 𝜇𝑘 ―𝑇 ∂𝑇 relation in Eq. (9) can be treated as a partial molar enthalpy, the other term within the square brackets is just a ratio of two kinetic coefficients, usually termed the heat of transport (Grout and Lidiard, 2008, Lidiard, 2015, Sugisaki, Sato and Furuya, 1981, Grout and Lidiard, 2008). This parameter has not yet been shown to have a thermodynamic representation or a uniformly accepted theoretical basis except for the case of an ideal gas. For historical reasons, the overall transport due to a thermal gradient, termed Soret diffusion (Rahman and Saghir, 2014), or thermodiffusion, is written as: 𝐽𝑠 = ― 𝐷𝑘𝑆𝑇𝑛𝑘∇𝑇 (10) where 𝐷𝑘 denotes the self-diffusion coefficient and 𝑆𝑇 denotes an experimentally determined Soret coefficient. Based on Eq. (9), the diffusion coefficient is embodied by the corresponding terms in the prefactor to the differential relations, whereas the Soret coefficient represents the complex expression within the square brackets. One of the difficulties in assigning a theoretical model for thermodiffusion is this mixed thermodynamic and kinetic nature. Recently developed diffusion formulations have proposed that thermodiffusion is entirely entropic in nature (Semenov and Schimpf, 2009, Kocherginsky, 2010). In addition, it has not been established how to obtain thermodiffusion and phenomenological proportionality coefficients from Eqs. (5) and (9) (Hartung and Kohler, 2009). 2.2 2.2.1 MECHANISTIC MODEL FOR MASS DIFFUSION Mass flux representation In the following we adopt an alternative approach for modeling mass diffusion that is based on a mechanistic description of the diffusion process (Kocherginsky and Gruebele, 2016, Eliaz and Banks5 Silis, 2008). Using the volumetric density of a component, 𝑛𝑘, as the primary variable with units of [mol m3], its flux is defined as the amount of component transported across a unit area normal to the flux direction, over a unit of time. Assuming that the diffusing component can be assigned an average constant drift velocity, 𝑣𝑘, its molar flux, 𝐽𝑘, is: 𝐽𝑘 = 𝑣𝑘 𝑛𝑘 (11) 𝑣𝑘 = 𝑀𝑘 𝐹𝑘 (12) where 𝐽𝑘 has units of mol m2s , and 𝑣𝑘 has units of m s. Using a linear approximation for uncorrelated movements of transporting particles, the drift velocity is a result of the product of the driving force, 𝐹𝑘, imparted to a particle, and the proportionality factor, 𝑀𝑘: The proportionality factor 𝑀𝑘, termed mechanical mobility, is an inverse of the frictional drag experienced by a particle as it interacts with its environment. The driving forces can be derived from thermal activation models for particle jump between neighboring sites. At the system scale, these forces are driving the system to the equilibrium state, which for constant temperature corresponds to the minimization of the Gibbs free energy, or in the case of nonuniform temperature, to the minimization of the Planck potential. The macroscopic model for the driving force towards uniform thermodynamic potential can therefore be written as (Einstein, 1926): which results in the mass flux equation: 𝐹𝑘 = ―∇𝜇𝑘 𝐽𝑘 = 𝑀𝑘𝑛𝑘𝐹𝑘 = ― 𝑀𝑘𝑛𝑘∇𝜇𝑘 (13) (14) The assumption that the velocity of the diffusing species is proportional to the gradient of chemical potential is also central to the transport model of irreversible thermodynamics. For simultaneous variation in chemical composition and temperature, recent models have proposed casting mass diffusion using a generalized physicochemical potential, 𝜇𝑔𝑘 (Kocherginsky and Gruebele, 2016). A thermodynamic gradient of 𝜇𝑔𝑘 thus includes all factors, including terms related to thermal and pressure gradients. Similarly, in (Semenov and Schimpf, 2009), the gradient in chemical potential is complemented by a pressure gradient term within a nonequilibrium thermodynamic approach. Andersson and Agren (Andersson and Agren, 1992) developed a formulation based on gradients of chemical potentials utilizing thermodynamic models in a CALPHAD (Computer Coupling of Phase Diagrams and Thermochemistry) (Spencer, 2008) approach to represent mobility. In our model, we use the mobility model for nonstoichiometric uranium dioxide from (Moore, Gueneau and Crocombette, 2013). 2.2.2 Conservation equation for mass transport The conservation equation for implementation in a FEM solver can be written: 6 ∂𝑛𝑘 ∂𝑡 + ∇ ∙ ( ― 𝑀𝑘𝑛𝑘∇𝜇𝑔𝑘) = 0 (15) The generalized version of chemical potential relations assume transport between infinitesimal volumes with uniform intensive variables. However, in the computational implementation, the values of the composition, temperature and chemical potential are interpolated over control volumes, and this interpolation allows us to work with general chemical potentials if they are evaluated within the same element as the intensive variables, as will be demonstrated in the discussion of FEM implementation of the transport equation. In UO2 LWR fuel, the two main types of oxygen transport are by interstitial (in hyper-stoichiometric fuel, UO2+x) and vacancy (in hypo-stoichiometric fuel, UO2-x) diffusion (Moore, Gueneau and Crocombette, 2013). In the case of interstitial diffusion, interstitial solute oxygen moves among interstitial oxygen sites in the host UO2 lattice. The defects that facilitate transport are controlled by stoichiometry, and, in the case of interstitial diffusion, the movement is not accompanied by a countermovement of a defect, such as a vacancy. If the underlying intrinsic mechanism of transport and the mobility for each thermodynamic driving force is the same, then such an interpolation can be further simplified by the use of a common mobility factor. In the case where the concentration of the moving species is given as a fraction of its total amount, such as the mole fraction of interstitial oxygen, 𝑦𝑖𝑂, related to the interstitial site fraction of, 𝑦𝑖, in the UO2+x, the flux in Eq. (14) has to be accordingly scaled by the same factor. 3. FEM IMPLEMENTATION OF THE MECHANICSTIC TRANSPORT EQUATIONS It is assumed that the temperature distribution is known and is treated as an externally imposed field, and thus the energy transport and balance equations will not be described here. This assumption is reasonable for nuclear fuel where the direct effect of mass transport on energy transport are negligible compared to the energy generated by fission, and influence heat conduction only by altering material thermal properties. Engineering transport models are based on spatial and temporal discretization of the domain of interest. They are usually cast in the form of control volumes and the integral form of conservation equations. The fuel performance code Bison uses a Lagrangian FEM approach in which the weak form of conservation equations is solved. The FEM approach tracks the material volume, which is more suitable for mechanical and thermomechanical analyses that are the focus of the fuel performance model, than species transport. To address species transport using the relations of irreversible thermodynamics described above, the discretized form of Eq. (15) uses mass concentration of the diffusing specie as the primary variable. It is based on local equilibrium within each infinitesimal element so that it can be described by local intensive variables, composition, and temperature. The FEM discretized form, however, uses the intensive variables at the nodes of the FEM mesh, and their values are interpolated using shape functions in the finite elements. The thermochemical equilibrium and resulting chemical potentials can be calculated not only at the nodes, but also at the FEM integration points using interpolated values of composition and temperature. 7 The two forms of the discretized processes are shown in Figure 2 (Niven and Noack, 2014). Δμ J Δ dz J dz dy μ dy dx μ2 μ1 dx (a) (b) Figure 2. Volume elements for (a) different, adjacent local equilibrium, and for (b) a continuum representation. The discretization in Figure 2a represents processes where each localized volume is at equilibrium for the given composition and temperature. The chemical potentials are thus at a minimum in each volume but individual species may have different values. Even though each volume is at local equilibrium, a physical flux occurs to diminish the spatial difference in chemical potentials. In this form, entropy production cannot occur in a volume but only between the volumes, which requires boundary entropy production terms and partition of entropy between the two volumes (Niven and Noack, 2014). In the second form, depicted in Figure 2b, each volume element does not have to be at equilibrium. And the composition and temperature can be represented as continuous variables at the nodes. Thermochemical equilibrium and the corresponding chemical potentials can be calculated at the nodes (vertices) of the volume, so that the gradients of chemical potentials are available within the volume. This form provides continuity of intensive variables and thermodynamic functions and is commonly adopted in fluid mechanics and heat transfer. However, it contradicts the assumption of equilibrium in a local volume, creating a philosophical difficulty in the use of intensive variables which are strictly defined only at equilibrium. When the equilibrium is enforced at the nodes, we can assume that the interpolated values of chemical potential within a volume are a good representation of the real values. The continuity of intensive variables also eliminates a need for entropy production at the volume boundaries (Niven and Noack, 2014). We deem the second form for implementation of equations Eqs. (14) and (15) more applicable to the problem at hand. The composition and temperatures are defined at the nodes where the thermodynamic equilibrium and the chemical potentials are also evaluated. Consequently, the gradient of chemical potential evaluated at the integration points of the finite elements corresponds to the total gradient as proposed in (Kocherginsky and Gruebele, 2016). The integral form of the mass conservation equation for implementation in the FEM framework is: 8 ∫ 𝑉 ∂𝑛𝑘 𝜓 𝑑𝑉 ― ∂𝑡 ∫𝐽 𝑉 𝑘 ∇𝜓 𝑑𝑉 + ∫𝜓 𝐽 Γ 𝑘 ∙ 𝑛 𝑑Γ = 0 (16) where V and Γ denote the volume and surface of the domain, 𝜓 denotes the test functions, and 𝑛 denotes the outward normal to the domain surface. In more compact notation, as used in the MOOSE (Gaston, Permann, Peterson, Slaughter, Andes, Wang, Short, Perez, Tonks, Ortensi, Zou and Martineau, 2015) documentation, the integral Eq. (16) is written as: 𝑅(𝑐𝑘) = ∂𝑛𝑘 ( ∂𝑡 ) , 𝜓 ― (∇𝜓 , 𝐽𝑘) + 〈𝜓 , 𝐽𝑘 ∙ 𝑛〉 = 0 (17) where parentheses and angled brackets denote corresponding inner products integrated across the domain and over its surface, respectively. In FEM, the concentration field, 𝑛𝑘, at any location is represented by values at the nodes of a finite element mesh, 𝑛𝑘𝑗, and their nodal shape functions, 𝜙𝑗, such that: 𝑛𝑘 = ∑𝑛 𝑘𝑗𝜙𝑗 𝑗 (18) In the FEM formulation used in MOOSE, the test and shape functions are the same. Using Eq. (18), the residual of Eq. (17) becomes a vector, 𝑹(𝒏𝒌), of the vector of nodal values, 𝒏𝒌, and the resulting system of nonlinear equations can be solved using Newton’s method by iteratively driving the residual vector to zero. The components of the residual vector are: 𝑅𝑖(𝒏𝒌) = ∂𝒏𝒌 ( ∂𝑡 ) , 𝜓𝑖 ― (∇𝜓𝑖 , 𝐽𝑘) + 〈𝜓𝑖 , 𝐽𝑘 ∙ 𝑛〉 , 𝑖 = 1,..,𝑁 (19) where index i denotes the finite element mesh node index associated with the test function 𝜓𝑖 , and N denotes the number of FEM nodes. The elements of the Jacobian matrix, 𝒥𝑖𝑗, required for the Newton iterations are: 𝒥𝑖𝑗(𝒏𝒌) = ∂𝑅𝑖 (𝒏𝒌) (20) ∂ 𝑛𝑘𝑗 where indices i and j denote nodes in the FEM mesh, and 𝑛𝑘𝑗 denotes a value of nk at node j. The inner products in Eq. (19) involve derivatives of the mass flux with respect to concentration at node j: ( ∇𝜓𝑖 , ∂𝐽𝑘 ∂𝑛𝑘𝑗 ) 〈 9 , 𝜓𝑖 , ∂𝐽𝑘 ∂𝑛𝑘𝑗 〉 ∙ 𝑛 (21) and refer to derivatives of the mass flux in the equation for 𝐽𝑘, which are not trivial. MOOSE’s Jacobianfree Newton Krylov solver emphasizes the importance of these calculations. Approximate expression for the Jacobian is, nevertheless, necessary in order to achieve a reasonable convergence rate. Using ∂𝒏𝒌 ∂𝑛𝑘𝑗 = ∂ ∑∂𝑛 𝑖 𝑘𝑗 (𝑛𝑘𝑖𝜙𝑖) = 𝜙𝑗 (22) for the mass flux in Eq. (17), the Jacobian inner product terms can be approximated as: (∇𝜓𝑖 , ― 𝑀𝑘𝐹𝑘𝜙𝑗) , 〈𝜓𝑖 , ― 𝑀𝑘𝐹𝑘𝜙𝑗 ∙ 𝑛〉 (23) which omits the partial derivatives of ― 𝑀𝑘 and 𝐹𝑘 . Eq. (14) is linear in the concentration variable, and the resulting conservation equation has advective form. This type of equation is notoriously difficult to solve and requires special time and space integration algorithms (Ewing and Wang, 2001). Recently, an upwind time integration scheme has been implemented and will be used for modeling the above formulated mass transfer using the calculated particle average drift velocity. 4. MODELING OF THE FUEL COMPOSITION WITH BURNUP Nuclear fuel composition evolves during burnup as actinides are consumed during fission and new elements are created through fission and transmutation. Thus, it is important to account for the effect of changing elemental composition on oxygen transport. Details of nuclear fuel burnup described in (Walker, Rondinella, Papaioannou, Van Winckel, Goll and Manzel, 2005) and modeled in (Piro, Banfield, Clarno, Simunovic, Besmann, Lewis and Thompson, 2013) were used for an example. The elemental composition along the radius of the pellet as a function of time was calculated using the Origen depletion code (Gauld, Hermann and Westfall, 2005). The composition was read into Bison as a Solution Function. It was assumed that oxygen is the only transporting element, and that it does not participate in nuclear reactions. Therefore, its total amount in the pellet remains constant, but its distribution may change due to mass diffusion. 5. IMPLEMENTATION OF TRANSPORT MODEL IN BISON The implementation of the oxygen transport formulation in Bison utilizes Thermochimica subroutines, Multiphysics Object Oriented Simulation Environment (MOOSE) functions, and MOOSE data structures. The Thermochimica subroutines are used for calculation of thermodynamic parameters that are used for the transport model and driving forces. The MOOSE FEM implementation is based on the weak form of conservation equations, namely mass flux and heat diffusion, and on constitutive material models and their parameters. The U-O system is used for simplified illustration. More realistic, complex systems can be modeled by using a larger elemental inventory as produced during burnup, with the dioxide treated as 10 MO2+x where fission product and transuranic elements substitute for uranium. MOOSE functions implemented for this purpose are: OxygenChemTransportAux (MOOSE AuxKernel) This function calculates the thermodynamic parameters required for oxygen transport properties. It uses temperature, pressure and elemental composition at the nodes of the FEM mesh, and calculates thermochemical equilibrium. It provides chemical potential of oxygen, Planck potential for the determined oxygen-to-uranium ratio, departure from stoichiometry in 𝑈𝑂2 ± 𝑥 , site fraction of oxygen vacancies on the regular anion UO2 sublattice, site fraction of oxygen interstitials on the oxygen interstitials sublattice (Moore, Gueneau and Crocombette, 2013), and the ratio of the oxygen interstitials to the total oxygen content. This function couples the thermodynamic models to parameters used in the heat and mass transport models. MobilityUOX (MOOSE Material) This function calculates the mobility parameter for uranium oxide based on (Moore, Gueneau and Crocombette, 2013). It uses site fractions of oxygen vacancies and interstitials as calculated in OxygenChemistryAux. It also couples the temperature field from the heat diffusion calculation. ChemicalDiffusion (MOOSE Kernel) This kernel uses the gradient in oxygen chemical potential from OxygenChemTransportAux as the driving force for the oxygen flux. The oxygen chemical potential gradient is multiplied by the oxygen mobility calculated in MobilityUOX, and by the concentration of oxygen multiplied by the ratio of the concentration of oxygen in the interstitial lattice to the total oxygen concentration. The main variable for this kernel is the total concentration of oxygen. OxygenDiffusionMobility (MOOSE Kernel) This kernel uses the gradient in oxygen concentration as the driving force for the oxygen flux. The oxygen concentration gradient is multiplied by the oxygen mobility calculated in MobilityUOX, and by the concentration of oxygen multiplied by the ratio of the concentration of oxygen in the interstitial lattice and the total oxygen concentration. The main variable for this kernel is the total concentration of oxygen. The relation between the BISON functions used for the transport simulation is shown in Table 1. For each function, the list of input (In) and updated (Out) variables determines the function’s dependency on other functions. Table 1. Relation between functions and variables in the oxygen transport simulation. Function MobilityUO2 OxygenChemTransportAux Variables 𝑦𝑉𝑎 , 𝑦𝐼𝑜 , T In Out MO In EL , O, T Out 𝑦𝑉𝑎 , 𝑦𝐼𝑜 , x , 𝜇𝑂 , 𝜇𝑂/𝑇, OIo/O Out OIo/O 11 OxygenDiffusionMobility In Out In ChemicalDiffusion Out In TimeDerivative Out O , T, OIo/O O O , 𝜇𝑂 , MO , OIo/O O O O Symbols in Table A1: Mo – mobility of oxygen, EL – Elemental composition, O – total concentration of oxygen atoms, T – temperature, 𝑦𝑉𝑎 - site fraction of oxygen vacancies on the regular anion sublattice, 𝑦𝐼𝑜 - site fraction of oxygen interstitial atoms, 𝜇𝑂 - chemical potential of oxygen, 𝜇𝑂/𝑇 – oxygen Planck potential, OIo – concentration of oxygen interstitial atoms, and x – departure from stoichiometry. Primary variables are denoted by boldface. 6. EXAMPLES Two example problems are presented. The first example is a model of oxygen transport in a prescribed temperature field, and composition history which results in a significant oxygen potential gradient, and resulting considerable oxygen redistribution. The second example applies the same conditions but uses the concentration gradient as the sole driving force for transport, and thus shows almost no oxygen redistribution. 6.1 OXYGEN TRANSPORT DRIVEN BY CHEMICAL POTENTIAL GRADIENT An axisymmetric model of a radial section of a fuel pellet from (Walker, Rondinella, Papaioannou, Van Winckel, Goll and Manzel, 2005) was used for the simulation of oxygen transport during burnup. Pellet diameter was 9.3mm and a parabolic temperature profile was imposed with a fuel centerline temperature of 1873.15K and fuel pellet surface temperature 673.15K. The material composition was calculated in (Piro, Banfield, Clarno, Simunovic, Besmann, Lewis and Thompson, 2013) and imposed as a function of time and radial location. Figure 3 shows the total power history and calculated final composition for the pellet. We ran our simulation only for the first 24 days as the purpose was to demonstrate the capability of the formulation. The equivalent composition at that time is shown in Figure 4. 12 Figure 3. Simulated elemental composition of the pellet at the end of the burnup. Figure 4. Simulated elemental composition of the after 24 days. The mobility model for non-stoichiometric uranium oxide from Reference (Moore, Gueneau and Crocombette, 2013) was used. The actual values were multiplied by 100 in order to speed up the transport process. Figures 5 shows the FEM mesh with overlaid with the temperature profile. 13 Figure 5. Temperature distribution across the FEM mesh along the pellet. Temperature [K] Figures 6-11 show the temperature, oxygen distribution, oxygen to uranium ratio, oxygen chemical potential calculated from the UO2±x phase (Eq. (7) in (Piro, Welland and Stan, 2015)), site fraction of oxygen interstitials, and the site fraction of oxygen vacancies along the radius of the pellet after 24 days (2.074*106 sec) of burnup. 1900 1800 1700 1600 1500 1400 1300 1200 1100 1000 900 800 Temperature 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 6. Temperature distribution in the pellet along the radius. 0.675 24 days O mol fraction [wrt t=0] 0.674 0.673 0.672 0.671 0.67 0.669 0.668 0.667 0.666 0.665 0 0.5 1 1.5 2 2.5 Radius [mm] 14 3 3.5 4 4.5 Figure 7. Oxygen distribution after 24 days in the pellet along the radius (normalized to mol fraction at the start time). 2.03 0 days 24 days 2.025 O/U 2.02 2.015 2.01 2.005 2 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 8. Oxygen to uranium ratio in the pellet along the radius at the start, and after 24 days. Oxygen Potential [kJ/mol] -240 0 days 24 days -260 -280 -300 -320 -340 -360 -380 -400 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 9. Oxygen chemical potential in the pellet along the radius at the start, and after 24 days. 15 0.03 0 days 24 days Site fraction of Oi 0.025 0.02 0.015 0.01 0.005 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 10. Site fraction of oxygen interstitials at the start, and after 24 days. 5 0 days 24 days Site fraction of Ova * 10 6 4.5 4 3.5 3 2.5 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 11. Site fraction of oxygen vacancies at the start, and after 24 days. Simulations show movement of oxygen down the gradient of the chemical potential towards the center of the pellet. 16 6.2 OXYGEN TRANSPORT DRIVEN BY CONCENTRATION GRADIENT In this example, we use concentration gradient as the driving force. The diffusion coefficient, 𝐷𝑘, is calculated from the Einstein relation (Einstein, 1926) 𝐷𝑘 = 𝑅𝑇𝑀𝑘, so that the driving flux is: 𝐽𝑘 = ― 𝑅𝑇𝑀𝑘∇𝑛𝑘 (24) 𝜇𝑘 = 𝜇𝑜 + 𝑅𝑇ln 𝑛𝑘 (25) which can be derived from equation (14) using chemical potential: Figure 12 shows the oxygen distribution after 24 days of burnup. Almost no transfer of oxygen occurred using the concentration model and the values in the graph are essentially numerical noise. The model would require addition thermal diffusion term from Eq. (10) to force the redistribution of the oxygen. O mol fraction [wrt t=0] 0.66726 0 days 24 days 0.66725 0.66725 0.66724 0.66724 0.66723 0.66723 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 13. Oxygen distribution in the pellet along the radius at the start, and after 24 days. Figures 14-15 show the oxygen to uranium ratio, and site fraction of oxygen interstitials along the radius of the pellet after 24 days (2.074*106 sec) of burnup. 17 2.0085 2.008 O/U 2.0075 2.007 0 days 24 days 2.0065 2.006 2.0055 2.005 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 14. Oxygen to uranium ratio in the pellet along the radius at the start, and after 24 days. 0.0085 Site fraction of Oi 0.008 0.0075 0.007 0 days 24 days 0.0065 0.006 0.0055 0.005 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Radius [mm] Figure 15. Site fraction of oxygen interstitials at the start, and after 24 days. The model for oxygen redistribution using Soret formalism has been described in a previous report and will not be repeated here. 18 7. SUMMARY We have developed a formulation for oxygen transport in urania nuclear fuels with compositional and temperature non-uniformity that more accurately represents energetic driven redistribution. The proposed model for oxygen transport is based on the full gradient of the oxygen chemical potential that is evaluated at the nodal points of the FEM mesh of the transport model. The formulation was demonstrated on the problem of temporally and spatially varying radial composition of nuclear pellet undergoing burnup. The computational implementation requires further development due to the advection character of the transport equation. Optimal strategy for defining material composition and improving the efficiency of thermodynamics calculations present additional opportunities for improvement of the modeling strategy. 8. REFERENCES Olander, D. R. (1976). Fundamental Aspects of Nuclear Reactor Fuel Elements, Technical Information Center, Office of Public Affairs, Energy Research and Development Administration. Olander, D. R. (1998). "Mechanistic interpretations of UO2 oxidation." Journal of Nuclear Materials, 252(1-2), 121-130. Walker, C. T., Rondinella, V. V., Papaioannou, D., Van Winckel, S., Goll, W., and Manzel, R. (2005). "On the oxidation state Of UO2 nuclear fuel at a burn-up of around 100 MWd/kgHM." Journal of Nuclear Materials, 345(2-3), 192-205. Park, K., Yang, M. S., and Park, H. S. (1997). "The stoichiometry and the oxygen potential change of urania fuels during irradiation." Journal of Nuclear Materials, 247, 116-120. Marchant, D. D., and Bowen, H. K. (1975). "Oxygen Redistribution in UO2 Due to a Temperature Gradient." Mass Transport Phenomena in Ceramics, A. R. Cooper, and A. H. Heuer, eds., Springer, 97-109. Adamson, M. G., Aitken, E. A., Evans, S. K., and Davies, J. H. (1975). "Oxygen Redistribution and its Measurement in Irradiated Oxide Fuels." Thermodynamics of Nuclear Materials 1974, IAEA, Vienna, Austria, 59-71. Adamson, M. G., and Carney, R. F. A. (1974). "Mechanistic Study of Oxygen Thermal-Diffusion in Hyperstoichiometric Urania and Urania-Plutonia Solid-Solutions." Journal of Nuclear Materials, 54(1), 121-137. Evans, S. K., Aitken, E. A., and Craig, C. N. (1969). "Effect of a Temperature Gradient on Stoichiometry of Urania-Plutonia Fuel." Journal of Nuclear Materials, 30(1-2), 57-61. Marin, J. F., and Contamin, P. (1969). "Uranium and Oxygen Self-Diffusion in UO2." Journal of Nuclear Materials, 30(1-2), 16 - 25. 19 Janek, J., and Timm, H. (1998). "Thermal diffusion and Soret effect in (U,Me)O2+delta: the heat of transport of oxygen." Journal of Nuclear Materials, 255(2-3), 116-127. Lassmann, K. (1987). "The Oxired Model for Redistribution of Oxygen in Nonstoichiometric Uranium Plutonium Oxides." Journal of Nuclear Materials, 150(1), 10-16. Olander, D. R. (1972). "Oxygen Redistribution in UO2+x." Journal of Nuclear Materials, 44(1), 116120. Matzke, H. (1987). "Atomic Transport-Properties in UO2 and Mixed Oxides (U, PU)O2." Journal of the Chemical Society-Faraday Transactions II, 83, 1121 - 1142. Lindemer, T. B., and Besmann, T. M. (1985). "Chemical Thermodynamic Representation of UO2+x." Journal of Nuclear Materials, 130(FEB), 473-488. Besmann, T. M., McMurray, J. W., and Simunovic, S. (2016). "Application of thermochemical modeling to assessment/evaluation of nuclear fuel behavior." Calphad-Computer Coupling of Phase Diagrams and Thermochemistry, 55, 47-51. Gueneau, C., Baichi, M., Labroche, D., Chatillon, C., and Sundman, B. (2002). "Thermodynamic assessment of the uranium-oxygen system." Journal of Nuclear Materials, 304(2-3), 161-175. Gueneau, C., Dupin, N., Sundman, B., Martial, C., Dumas, J. C., Gosse, S., Chatain, S., De Bruycker, F., Manara, D., and Konings, R. J. M. (2011). "Thermodynamic modelling of advanced oxide and carbide nuclear fuels: Description of the U-Pu-O-C systems." Journal of Nuclear Materials, 419(13), 145-167. Besmann, T. M. (2012). "Computational Thermodynamics: Application to Nuclear Materials." Comprehensive Nuclear Materials, R. J. M. Konings, ed., Elsevier, 455 - 470. Demirel, Y. (2007). Nonequilibrium Thermodynamics, Elsevier. Tschoegl, N. W. (2000). Fundamentals of Equilibrium and Steady-State Thermodynamics, Elsevier Science & Technology Books. Piro, M. H. A., Simunovic, S., Besmann, T. M., Lewis, B. J., and Thompson, W. T. (2013). "The thermochemistry library Thermochimica." Computational Materials Science, 67, 266-272. Piro, M. H. A., Banfield, J., Clarno, K. T., Simunovic, S., Besmann, T. M., Lewis, B. J., and Thompson, W. T. (2013). "Coupled thermochemical, isotopic evolution and heat transfer simulations in highly irradiated UO2 nuclear fuel." Journal of Nuclear Materials, 441(1-3), 240-251. Matzke, H. (1981). "Diffusion in Nonstoichiometric Oxides." Nonstoichiometric Oxides, T. S. Sorensen, ed., Academic Press, 155 - 269. de Groot, S. R., and Mazur, P. (2011). Non-Equilibrium Thermodynamics, Dover Publications, Inc. Kocherginsky, N., and Gruebele, M. (2016). "Mechanical approach to chemical transport." Proceedings of the National Academy of Sciences of the United States of America, 113(40), 11116-11121. Moore, E., Gueneau, C., and Crocombette, J.-P. (2013). "Diffusion model of the non-stoichiometric uranium dioxide." Journal of Solid State Chemistry, 203, 145-153. Hales, J. D., Novascone, S. R., Spencer, B. W., Williamson, R. L., Pastore, G., and Perez, D. M. (2014). "Verification of the BISON fuel performance code." Annals of Nuclear Energy, 71, 81-90. 20 Zienkiewicz, O. C., Taylor, R. L., and Zhu, J. Z. (2013). The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann. Hillert, M. (2007). Phase Equilibria, Phase Diagrams and Phase Transformations, Their Thermodynamic Basis, Cambridge University Press, New York. Kondepudi, D., and Prigogine, I. (2015). Modern thermodynamics: from heat engines to dissipative structures, John Wiley & Sons, Ltd. Onsager, L. (1931). "Reciprocal Relations in Irreversible Processes I." Physical Review, 37, 405-426. Grout, P. J., and Lidiard, A. B. (2008). "Computation of heats of transport in crystalline solids: II." Journal of Physics-Condensed Matter, 20(42). Lidiard, A. B. (2015). "An Essay on the Heat of Transport in Solids and a Partial Guide to the Literature." Diffusion Phenomena in Engineering Materials, I. Belova, G. Murch, and A. Öchsner, eds., 57-68. Sugisaki, M., Sato, S., and Furuya, H. (1981). "Solid-Phase Thermal-Diffusion as a Mechanism for Oxygen Redistribution in UO2+/-x and (U,PU)O2+/-x." Journal of Nuclear Materials, 97(1-2), 7987. Grout, P. J., and Lidiard, A. B. (2008). "Computation of heats of transport in crystalline solids: II." Journal of Physics-Condensed Matter, 20(42), 1-7. Rahman, M. A., and Saghir, M. Z. (2014). "Thermodiffusion or Soret effect: Historical review." International Journal of Heat and Mass Transfer, 73, 693-705. Semenov, S. N., and Schimpf, M. E. (2009). "Mass transport thermodynamics in nonisothermal molecular liquid mixtures." Physics-Uspekhi, 52(10), 1045-1054. Kocherginsky, N. (2010). "Temperature-driven mass, charge and heat transfer in terms of physicochemical potential and Einstein's mobility." Chemical Engineering Science, 65(14), 41544159. Hartung, M., and Kohler, W. (2009). "Reversible mass exchange between two multicomponent systems of different temperatures." European Physical Journal E, 29(1), 117-121. Eliaz, N., and Banks-Silis, L. (2008). "Chemical potential, diffusion and stress - Common confusions in nomenclature and units." Corrosion Reviews, 26(2-3), 87-103. Einstein, A. (1926). Investigations on the theory of the Brownian motion, Dover Publications. Andersson, J. O., and Agren, J. (1992). "Models for Numerical Treatment of Multicomponent Diffusion in Simple Phases." Journal of Applied Physics, 72(4), 1350-1355. Spencer, P. J. (2008). "A brief history of CALPHAD." Calphad-Computer Coupling of Phase Diagrams and Thermochemistry, 32(1), 1-8. Niven, R. K., and Noack, B. R. (2014). "Control volume analysis, entropy balance and the entropy production in flow systems." Beyond the Second Law: Entropy Production and Non-Equilibrium Systems, R. C. Dewar, C. Lineweaver, R. K. Niven, and K. Regenauer-Lieb, eds., Springer-Verlag, 129-162. 21 Gaston, D. R., Permann, C. J., Peterson, J. W., Slaughter, A. E., Andes, D., Wang, Y., Short, M. P., Perez, D. M., Tonks, M. R., Ortensi, J., Zou, L., and Martineau, R. C. (2015). "Physics-based multiscale coupling for full core nuclear reactor simulation." Annals of Nuclear Energy, 84, 45-54. Ewing, R. E., and Wang, H. (2001). "A summary of numerical methods for time-dependent advectiondominated partial differential equations." Journal of Computational and Applied Mathematics, 128(1-2), 423-445. Gauld, I., Hermann, O., and Westfall, R. (2005). "ORIGEN Scale System Module to Calculate Fuel Depletion, Actinide Transmutation, Fission Product Buildup and Decay, and Associated Radiation Terms." Piro, M. H. A., Welland, M. J., and Stan, M. (2015). "On the interpretation of chemical potentials computed from equilibrium thermodynamic codes." Journal of Nuclear Materials, 464, 48-52. 22