An Eulerian Single-Phase Transport Model for Solid Fission Products in the Molten Salt Fast Reactor: Development of an Analytical Solution for Verification Purposes

Nuclear reactor modeling has been shifting, over the last decades, towards full-core multiphysics analysis due to the ever-increasing safety requirements and complexity of the designs of innovative systems. This is particularly true for liquid-fuel reactor concepts such as the Molten Salt Fast Reactor (MSFR), given their strong intrinsic coupling between thermal-hydraulics, neutronics and fuel chemistry. In the MSFR, fission products (FPs) are originated within the liquid fuel and are carried by the fuel flow all over the reactor core and through pumping and heat exchange systems. Some of FP species, in the form of solid precipitates, can represent a major design and safety challenge, e.g., due to deposition on solid boundaries, and their distribution in the core is relevant to the design and safety analysis of the reactor. In this regard it is essential, both for the design and the safety assessment of the reactor, the capability to model the transport of solid FPs and their deposition to the boundary (e.g., wall or heat exchanger structures). To this aim, in this study, models of transport of solid FPs in the MSFR are developed and verified. An Eulerian single-phase transport model is developed and integrated in a consolidated multiphysics model of the MSFR based on the open-source CFD library OpenFOAM. In particular, general mixed-type deposition boundary conditions are considered, to possibly describe different kinds of particle-wall interaction mechanisms. For verification purposes, analytical solutions for simple case studies are derived ad hoc based on the extension of the classic Graetz problem to linear decay, distributed source terms and mixed-type boundary conditions. The results show excellent agreement between the two models, and highlight the effects of decay and deposition phenomena of various intensity. The resulting approach constitutes a computationally efficient tool to extend the capabilities of CFD-based multiphysics MSFR calculations towards the simulation of solid fission products transport.


INTRODUCTION
Liquid-fuel reactor concepts have gained renewed interest over the last years. Among them, the Molten Salt Reactor (MSR) and, in particular, the fast-spectrum MSFR (Molten Salt Fast Reactor) has obtained a prominent role thanks to the selection as one of the Generation IV reference technologies (Serp et al., 2014). The adoption of a circulating liquid fuel, in conjunction with the fast neutron spectrum, makes the MSFR system unique from the design and modeling viewpoints. Internal heat generation, fuel thermal feedback and transport of delayed neutron precursors and fission products lead to a strong intrinsic coupling between thermal-hydraulics, neutronics and fuel chemistry. Reactor modeling efforts have therefore shifted towards full-core and multiphysics analysis to meet the requirements and complexity of physical and computational models for the MSFR. A comprehensive account of state-of-the-art multiphysics modeling tools for the MSFR can be found in the work of Tiberga et al. (2020).
In the MSFR, fission products (FPs) that are originated within the fuel as the result of fission reactions are not retained by solid fuel elements and are therefore free to be carried by the liquid fuel flow within the primary circuit. The presence of mobile FPs represents a major design and safety challenge, as some of them are not expected to form stable compounds with the constituents of the fuel salt mixture (Grimes, 1970;Baes, 1974) and therefore give rise to separate phases, either in the form of solid precipitates or gas bubbles. In the case of solid FPs, precipitates are likely to deposit on reactor solid surfaces (Kedl, 1972;Compere et al., 1975) with potential issues such as formation of localized decay heat sources as well as deterioration of heat exchanger performance. Surface deposits might also pose a serious radiological threat in inspection/maintenance operations. Despite the analysis of FPs transport in recent MSFR studies is still limited, information on their distribution over the reactor core has great relevance for the design of the reactor, and in particular for the analysis of removal/reprocessing (Delpech et al., 2009) and bubbling (Cervi et al., 2019b) units. In this regard, there has been focus on modelling the behavior and effects of Xenon transport in MSRs (Price et al., 2020). Concerning metallic FPs, their analysis in the context of a multiphysics system code for the study of the Molten Salt Reactor Experiment has been recently addressed by Walker and Ji (2021). The aim of this work is the development of models of transport of solid FPs and their integration in state-of-the-art multiphysics tools adopted for MSFR analysis. To this aim, a single-phase Eulerian transport modeling approach is chosen, which represents a common and versatile choice for particle transport modeling in complex turbulent flows, and allows for the direct implementation in consolidated MSFR codes. A similar approach has been recently employed for the coupling of CFD and equilibrium chemistry calculations with application to localized-corrosion modeling in lead-cooled reactors (Marino et al., 2020). The development platform is the C++ open-source finite-volume CFD library OpenFOAM (OpenFOAM, 2021).
A preliminary implementation is here described and verified against analytical solutions for simplified test cases. The analytical solutions are derived from the well-known Graetz problem, which is here formulated for a general case with distributed source terms, linear decay and mixed-type boundary conditions. To the authors knowledge, such form of the Graetz problem has not been addressed in detail in the literature, and therefore the complete derivation is here reported. The influence of some key model parameters is discussed, and corresponding solutions from the OpenFOAM and analytical models are presented and compared.
The paper is organized as follows. The adopted multiphysics approach is briefly described in Section 2.1. In Section 2.2, the proposed FPs transport model implemented in the OpenFOAM solver is described in detail. The analytical solution against which the OpenFOAM FPs transport model is verified is presented in Section 2.3. Section 3 presents the results of the verification together with some discussion. Conclusive remarks are finally reported in Section 4.

The Multiphysics Model
The OpenFOAM library, based on standard finite-volume methods for CFD calculations, is used to develop the numerical solver used in the present work. Originally developed for the transient analysis of the MSFR (Aufiero et al., 2014), the adopted multiphysics solver was recently extended to allow for the study of compressibility effects during super-prompt-critical transients (Cervi et al., 2019b) and of the bubbling system (Cervi et al., 2019a). The version employed in this work features single-phase incompressible thermal-hydraulics, multi-group neutron diffusion and transport equations for delayed neutron and decay heat precursors. Transport equations for fission products are solved alongside the other physical modules, to provide a fully-coupled multiphysics simulation. All model equations are solved by means of finite-volume discretization, following an iterative segregated coupling approach.

Thermal-Hydraulics Model
Continuity, momentum and energy (in temperature form) conservation equations are expressed in a single-phase incompressible formulation: where the quantities u, p and T, which correspond to velocity, pressure and temperature, respectively, are intended as averaged in the sense of Reynolds-Averaged Navier-Stokes modeling. It follows that turbulence modeling is performed by means of standard eddy-viscosity based closure models, such as the well-known k-ε model (Launder and Spalding, 1974), for which effective momentum and thermal diffusivities can be expressed as the sum of a laminar and a turbulent contribution: where Pr and Pr t are the Prandtl and turbulent Prandtl numbers, respectively. Momentum and energy equations are coupled thanks to the Boussinesq approximation, for which the density value driving the buoyancy term in Eq 2 is linearized around a reference temperature T 0 and β T represents the volumetric thermal expansion coefficient of the fluid. Except for what concerns density in the linearized buoyancy term, constant average values are used for thermophysical properties to keep the numerics and the coupling between different physics as simple as possible. Finally, q represents a volumetric energy source which includes internal heat generation (both prompt and delayed) and optionally other energy sinks to model heat removal systems. Pressure-velocity coupling is performed through the standard SIMPLE/PISO algorithms (Patankar and Spalding, 1972;Issa, 1986).

Neutronics Model
The multi-group diffusion model is adopted for neutron flux calculations (Hébert, 2010). Despite some limitations, it is widely employed in standard nuclear reactor analysis. Thanks to its relative simplicity and limited computational effort, it has found several successful applications especially for multiphysics analysis (Aufiero et al., 2014;Fiorina et al., 2016). More recent works have also proposed the extension of OpenFOAM-based multiphysics codes to more advanced neutron transport approaches, e.g. the SP3 model (Fiorina et al., 2017;Cervi et al., 2019b). The diffusion equation for the gth group-integrated neutron flux φ g reads: where S n,g is the explicit neutron source of the gth group, constituted by prompt fission and scattering neutrons from other groups and delayed neutron precursors decay: The presence of φ h in the explicit source term S n,g couples the transport equations for different energy groups, which are therefore dealt with following a segregated iterative approach.
Most symbols have straightforward meaning and are listed in the Nomenclature section. It is worth mentioning that k eff acts as a tunable multiplication factor to model a prescribed reactivity insertion. A power-iteration routine based on the k-eigenvalue method is also included for steady-state simulation, which allows for the iterative adjustment of k eff to attain criticality at a specified power level. Due to the circulating nature of the fuel, transport equations are formulated also for delayed neutron and decay heat precursors. The transport equation for the concentration of delayed neutron precursors of the kth family c k reads: A discussion on the diffusion coefficient D eff is given in Section 2.2. We omit here for brevity the transport equation for the decay heat precursors, which is entirely analogous.
Group constants are adjusted as functions of local temperature around reference values to account for Doppler and fuel density effects. For a generic neutron reaction r occurring in the gth energy group: where Doppler effects are modeled by means of a logarithmic term where Σ 0 r,g and A 0 r,g , respectively represent the cross-section and a corresponding logarithmic coefficient at a reference temperature T Σ 0 , whereas density effects are taken into account through a linear correction consistently with the buoyancy term. The reference temperature for cross-sections can be chosen independently from T 0 . An analogous approach is employed for the correction of the intra-group neutron diffusion coefficient D n,g . The quantities Σ 0 r,g and A 0 r,g are evaluated by means of the Monte Carlo reactor physics and burnup code SERPENT 2 (Leppänen et al., 2015).

Fission Products Transport Model
Similarly to other transported scalar quantities, each fission product specie is modelled as a continuous scalar concentration field subject to advection, dispersion and decay mechanisms: where c is the concentration of the species under consideration, expressed in number of particles per unit volume. D eff is the total particle diffusivity, u is the carrier velocity field and λ and S represent the decay constant and source of a fission product specie, respectively. When chemical interactions between species and formation of separate phases are neglected, the source term can be simply related to the fission rate through a suitable yield coefficient y f : As previously mentioned, this modeling choice is motivated by the need to limit the overall complexity and computational requirements of the MSFR, and to easily integrate such models in state-of-the-art MSFR codes. The single-phase Eulerian approach can still represent a valid approximation, provided that some conditions are met. Theoretical and experimental analysis has suggested that Fick's diffusion law only applies when inertial effects are negligible, and that particles inertia plays an increasingly dominant role in transport mechanisms as particles size increases (Liu and Agarwal, 1974;Guha, 2008). Little information is known about the expected size of fission product particles in the MSFR, but previous experience with MSRs suggests formation of colloidal suspensions (i.e., with particle diameters approx. between 10 -9 and 10 -6 m) should be expected (Compere et al., 1975) and that simple diffusion laws may apply, at least as a first approximation. The other main condition requires the particle concentration in the carrier fluid to be sufficiently low to make all interactions between each particle and the fluid, or among particles themselves, negligible. Provided that such conditions are met, the adoption of more sophisticated approaches, e.g., Eulerian-Eulerian multiphase or Eulerian-Lagrangian, may prove little benefit at the expense of a much more complex modeling framework.
As regards the particle diffusivity, it is commonly assumed that the diffusivity coefficient D eff may be separated in a laminar and a turbulent contribution, analogously to Eq 5: where Sc and Sc t are Schmidt and turbulent Schmidt numbers, respectively. An alternative expression for the laminar diffusivity D is given by the so-called Einstein equation: which is derived under the assumption of large Schmidt number Sc (Balboa Usabiaga et al., 2013). k B is the Boltzmann constant and d p is the particle diameter. For monodisperse particles in isothermal bulk flow, D can therefore be regarded as a constant. As already stated, the inverse proportionality from Eq 13 between the particle flux and d p only holds for small values of d p , even in full-laminar flows.

Deposition Modeling
Besides particle transport in the bulk flow, transport mechanisms which lead to deposition need to be addressed separately. First of all, when particle-wall interaction in the boundary layer is considered, a variable diffusion coefficient can be introduced to model hydrodynamic interactions between particles and solid walls (Brenner, 1961). In such case, D eff is assumed as a function of the particle-wall distance: where y denotes the wall-distance in the normal direction in a boundary layer flow. Moreover, it is commonly accepted that D eff is mostly constant in the bulk of the flow and abruptly decreases in a very small layer close to walls. For this reason it can be assumed where D ∞ eff is a constant bulk diffusivity and F D (y) a correction factor which is always comprised within 0 and 1 (Brenner, 1961). Moreover, to model deposition mechanisms and formulate appropriate boundary conditions for the particle concentration field, one possible approach is the inclusion in the transport equation of a particle-wall interaction forcing term based on an interaction potential energy ϕ (Ruckenstein and Prieve, 1973;Bowen et al., 1976). Since boundary layer flows are inherently two-dimensional, close to a generic wall the steady-state transport equation reads where x and y here denote the longitudinal and transversal directions. This approach is general but introduces significant complication in the general problem. It is therefore beneficial, when possible, to decouple the advection-diffusion and deposition problems. When the energy potential ϕ and the hydrodynamic interaction effects are significant only over distances of the order of the particle size, this can be effectively achieved by dividing the boundary layer in two regions (Figure 1): the bulk region, where advection occurs and interactions are negligible, and a very thin wall region, where advection is negligible and deposition processes take place. The two regions are subject to an interface condition, which equates the particle concentration c(δ ϕ ) and its flux at a certain distance δ ϕ from the wall (with δ ϕ chosen such that ϕ(δ ϕ ) is arbitrarily small). It has been shown that, under fairly general hypotheses on ϕ, the coupled solution of the two regions may be approximated such that the wall region influence is collapsed in a first-order reaction boundary condition for the bulk region (Ruckenstein and Prieve, 1973;Prieve and Ruckenstein, 1976): The coefficient c is therefore a constant depending on F and ϕ. In principle its value might also depend on the integration limit δ ϕ , but it has been shown that in many circumstances its influence on the value of the integral is relatively small over a relatively broad range (Prieve and Ruckenstein, 1976). It should be noted that these results are derived under the assumption of negligible axial diffusion (which is commonly the case) and in absence of source terms and decay. They are, however, representative of a broad class of particle-wall interaction problems and it is here assumed that internal sources or decay phenomena do not alter significantly such behavior.
In the general three-dimensional case, the wall boundary condition can be written as for any point on the wall boundary, where n denotes the outward pointing wall-normal direction. Since particle-wall interactions are modeled with a simple first-order boundary condition, the deposited quantities can easily be tracked by solving where c d represents the surface concentration of deposited particles on the wall boundary, expressed in number of particles per unit area. In this simple model, wall adsorption of FP particles is modeled through a single deposition parameter c, which has the physical dimensions of a velocity. Desorption mechanisms can be included as well by adding a corresponding term Eq 20, with δ being the desorption rate constant (Wood et al., 2004).

Test Case
In this Section the results of the verification of the implemented FP transport model are described. A simple test case with an analytical solution has been chosen. A two-dimensional channel between parallel plates is considered (Figure 2). The particle transport model is decoupled from neutronics and heat transfer models. The source term S from Eq 10 is therefore specified explicitly, and a steady-state isothermal laminar flow is considered. The problem here considered resembles the well-known Graetz problem, for which different solutions are available in the literature. An exhaustive treatment of the Graetz theory applied to particle transport problems is given by Bowen et al. (1976), although it doesn't consider linear decay or distributed source terms. It is indeed difficult to find, in the literature, realistic applications which consider simultaneously, as in the case under consideration, distributed internal sources, decay reactions and mixed type boundary conditions. Corresponding analytical solutions for this specific problem are not found in the literature, and are therefore derived in the following.

Momentum Equation
Analytical solutions of the momentum equation can be found only for simple steady-state fully-developed laminar flow problems. In such a case, the well known parabolic solution reads u u y 0 with FIGURE 1 | Decoupling of the transport and deposition problems: solutions for the wall region (c w ) and the bulk region (c w ) are obtained separately and joined by imposing interface continuity on particle concentration and particle flux at y δ ϕ . The wall region thickness δ ϕ ≪ H is defined as the distance at which ϕ(δ ϕ ) becomes arbitrarily small. where u m is the maximum profile velocity and H is half the channel width as depicted in Figure 2.

Particle Concentration Equation
The boundary value problem then becomes, in explicit cartesian coordinates, where D is used, from now on, to denote the constant bulk diffusivity. Boundary conditions for the x-direction are discussed later. The problem is conveniently re-scaled by defining appropriate non-dimensional quantities: where c 0 is a concentration value typical of the problem. The equations are then rewritten as where non-dimensional groups are defined as Longitudinal diffusion is neglected to allow for separation of variables. This assumption is reasonable in all cases where diffusion is negligible compared to advection, i.e., if Pe ≫ f . Since the order of the equation with respect to x is now reduced, a single (inlet) boundary condition for x is required. As a further simplification, full symmetry of the problem with respect the x-axis is assumed. The full problem now reads where c in specifies a uniform inlet condition for c(0, y).
For the following discussion, it is convenient to assume that the distributed source can be expressed as where S is a representative value of the source, e.g., its average value over the domain. The choice of c 0 is purely a matter of convenience. A meaningful definition, however, is not trivial to find in the general case, given the interplay of several physical phenomena. In the simplest case with no internal source, concentration profiles become self-similar far from the inlet but no fully-developed solutions can be attained (in presence of removal mechanisms such as deposition and decay). In such case, the definition is straightforward: On the other hand, when a distributed source is present, the inlet contribution is forgotten as the fully-developed concentration profile is attained and therefore a more meaningful choice should be based on the relative intensity of generation and removal mechanisms. When radioactive decay is dominant, a good definition reads When decay is negligible, the reference concentration c 0 can be chosen as These values are useful to identify correct scaling with respect to the dominant removal mechanisms. Similar expressions are easily found when solving for the centerline concentration in fully-developed profiles with uniform source.
Solutions of the boundary value problem Eq 26 can be found by separation of variables:

Derivation of the Analytical Solution
The functions Y n ( y) coincide with the eigenfunctions of the associated Sturm-Liouville problem found by isolating the y part of the equation. This can be highlighted by inserting Eq 31 in Eq 26 and after some manipulation: The eigenvalue problem is therefore stated as where and β 2 n being the eigenvalue associated to the eigenfunction Y n . The problem must be equipped with the corresponding boundary conditions: It is easily verified that β n 0 is not an eigenvalue. Then, through the simple change of variable η 2β n y, Eq 35 can be recast into one form of the parabolic cylinder equation (DLMF, 2020, Ch. 12): Y n y A 1,n Y 1,n y + A 2,n Y 2,n y where Y 1,n y e −β nŷ 2 / 2 M 1 − β n 4 + Da 4β n , 1 2 ; β n y 2 (38) Y 2,n y 2β n y e −β nŷ 2 / 2 M 3 − β n 4 + Da 4β n , 3 2 ; β n y 2 (39) are two linearly independent solutions. M(·, ·; ·) is the confluent hypergeometric function of the first kind (DLMF, 2020, Ch. 13), which can be expressed as the following power series expansion: where (a) k a(a + 1)/(a + k − 1) denotes the so-called Pochhammer symbol, or rising factorial. The eigenvalue problem must be solved by imposing the associated boundary conditions. By means of the following properties of M(a, b; z), it can be easily shown that the symmetry condition at y 0 implies A 2,n 0. The wall boundary condition requires which gives the β n values that correspond to the non-trivial solutions of the eigenvalue problem. Solutions of Eq 44 have to be found numerically; more detailed numerical considerations can be found in the Supplementary Appendix. The eigenfunctions are therefore determined (up to an arbitrary multiplicative constant) as Y n y Y 1,n y (45) whose general solution (back in terms of y) can be expressed as The separated Eq 32 now reads The eigenfunctions Y n are orthogonal with respect to the scalar product defined as This is exploited to transform Eq 46 into To obtain this last form, S has been expanded as and the arbitrary constants X n (0) have been determined from the remaining inlet boundary condition (expanded as well): where FIGURE 3 | Dependence of the 1st eigenvalue β 1 on model parameters Sh and Da.
The normalization constants C n and D n can be found by numerical integration of Eqs 49, 55, respectively. More details on their computation are reported in the Supplementary Appendix.

RESULTS AND DISCUSSION
In this Section, results of the verification of the implemented models against the analytical solutions are presented. In Section 2.3.3 it is shown that, if the effect of desorption is negligible compared to decay, the β n are the roots of   Sh Y n (1) + Y ′ n (1) with Y n y e −β nŷ 2 / 2 M 1 − β n 4 + Da 4β n , 1 2 ; β n y 2 (57) It is therefore evident that the inclusion of a linear decay term in the transport equation affects the concentration profiles shape by shifting the eigenvalues. As shown in Figures 3-6, the influence of the decay parameter Da is significant. For dominant modes, the increase of Da tends to flatten the β n curves, resulting in a vanishing influence of Sh in determining the shape of the concentration field. On the contrary, higher-order modes are less affected, with the eigenvalue shift due to linear decay decreasing as n increases. In the following, results from the comparison between the OpenFOAM model described in Section 2.2 and the corresponding analytical solutions are discussed. To highlight the role of decay and deposition phenomena, the selected parameters are Da and Sh, with values ranging in [0.1,1,10] for both. Case study parameter values are listed in Table 1, while other relevant parameters of the problem are listed in Table 2.
To simplify the analysis, the inlet concentration is set to zero (c in 0). Furthermore, to allow for a direct comparison between different test cases, the reference concentration value has been selected as It follows that S x, y S L x, H y S S L x, H y Therefore in dimensionless form the solution does not depend on the source term average value, but only on its shape. For the present analysis, a cosine-shaped source term has been selected to resemble the typical shape of fission rate profiles in simplified reactor geometries such as the one here considered: Concentration profiles obtained for the nine test cases ( Table 1) are shown in Figures 7-9.
Results show excellent agreement between the proposed transport model and the analytical solutions, proving a successful verification of the implemented transport models in OpenFOAM. The influence of decay is evidenced from the decrease in concentration profiles from Da 0.1 to Da 10. Besides the average value, the effect of decay is also evident on the shape itself as also pointed out by the analysis of the eigenvalues. With increasing Da, particles can diffuse less before decaying, and therefore the concentration profiles tend to resemble more the shape of the source term. The same effect can be seen on the combined effect of deposition: as previously said, the effect of Sh on the profiles is more evident for smaller values of Da, while profiles tend to become more similar as Da increases.
We finally report here some brief computational information regarding the verification. All simulations were performed on a structured orthogonal mesh, constituted by 5 × 10 4 hexahedral volumes, with respectively 500 and 100 divisions on the longitudinal and transversal directions. All 9 cases shown similar convergence behavior, with approximately 5 × 10 4 pseudo-transient iterations needed to ensure tight convergence for the particle concentration field. Computational times are comparable among all cases, showing no significant dependence on the physical parameters within the selected range (Table 3).

CONCLUSION
In this paper, the preliminary development of transport models for the solid fission products in the MSFR was discussed. Simplified transport models based on a Eulerian single-phase framework were implemented in consolidated MSFR multiphysics simulation tools based on the open-source finite-volume library OpenFOAM. The resulting model has been verified against analytical solutions for simplified test cases, based on an extended version of the Graetz problem. Results show the good agreement between the two models, proving the correct implementation of the transport and deposition mechanisms considered and the capability of OpenFOAM in treating coupled deposition and decay phenomena of different relative intensities, albeit on a simplified reference problem. The proposed approach constitutes a computationally efficient framework to extend the capabilities of CFD-based multiphysics MSFR calculations towards the simulation of solid fission products transport.
The analytical model used for verification has been developed specifically for this application. The simultaneous presence of distributed internal generation, radioactive decay and mixed deposition boundary condition in a Graetz problem represents an original contribution, and the resulting analytical model could therefore constitute a useful benchmark for future developments of the FPs migration model and for other similar MSFR applications. Future work will include the analysis of the integration of FPs transport in a full reactor simulation. In this regard, turbulent transport in more complex geometries may lead to large concentration gradients close to walls, with the need for mesh refinement. The extension to reactor simulation might therefore lead to numerical and/or computational issues to be addressed. Among other possible limitations of the current modelling approach, we highlight the adoption pure concentration-driven diffusive transport. As already mentioned, such approximation is strictly valid only for particles of very small size. Further study will be needed to assess the validity of such assumption for MSFR applications and to possibly extend the methodology to more advanced particle transport models.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material. The data that support the findings of this study are openly available in Zenodo at http://doi.org/10. 5281/zenodo.4557568.