Skip to main content

ORIGINAL RESEARCH article

Front. Energy Res., 29 June 2021
Sec. Nuclear Energy
This article is part of the Research Topic Reactor Fuels, Materials and Systems under Extreme Environments View all 29 articles

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

  • 1Department of Energy, Politecnico di Milano, Milano, Italy
  • 2Sezione di Milano-Bicocca, Istituto Nazionale di Fisica Nucleare (INFN), Milano, Italy

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.

1 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.

2 Materials and Methods

2.1 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.

2.1.1 Thermal-Hydraulics Model

Continuity, momentum and energy (in temperature form) conservation equations are expressed in a single-phase incompressible formulation:

u=0(1)
ut+(uuT)=1ρp+[1βT(TT0)]g+[νeff(u+(u)T)](2)
Tt+(uT)=(αeffT)+qρcp(3)

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:

νeff=ν+νt(4)
αeff=α+αt=νPr+νtPrt(5)

where Pr and Prt 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 T0 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).

2.1.2 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:

1vgφgt=(Dn,gφg)Σa,gφghgΣs,ghφg+(1βd)χp,gν¯gkeffΣf,gφg+Sn,g(6)

where Sn,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:

Sn,g=(1βd)hgχp,hν¯hkeffΣf,hφh+hgΣs,hgφh+χd,hkλkck(7)

The presence of φh in the explicit source term Sn,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 keff 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 keff 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 ck reads:

ckt+(uck)=(Deffck)λkck+βd,kgν¯gΣf,gφg(8)

A discussion on the diffusion coefficient Deff 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:

Σr,g=(Σr,g0+Ar,g0logTT0Σ)1βT(TT0)1βT(T0ΣT0)(9)

where Doppler effects are modeled by means of a logarithmic term where Σr,g0 and Ar,g0, respectively represent the cross-section and a corresponding logarithmic coefficient at a reference temperature T0Σ, 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 T0. An analogous approach is employed for the correction of the intra-group neutron diffusion coefficient Dn,g. The quantities Σr,g0 and Ar,g0 are evaluated by means of the Monte Carlo reactor physics and burnup code SERPENT 2 (Leppänen et al., 2015).

2.2 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:

ct+(uc)=(Deffc)λc+S(10)

where c is the concentration of the species under consideration, expressed in number of particles per unit volume. Deff 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 yf:

S=yfgΣf,gφg(11)

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 Deff may be separated in a laminar and a turbulent contribution, analogously to Eq 5:

Deff=D+Dt=νSc+νtSct(12)

where Sc and Sct are Schmidt and turbulent Schmidt numbers, respectively. An alternative expression for the laminar diffusivity D is given by the so-called Einstein equation:

D=kBT3πρνdp(13)

which is derived under the assumption of large Schmidt number Sc (Balboa Usabiaga et al., 2013). kB is the Boltzmann constant and dp 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 dp only holds for small values of dp, even in full-laminar flows.

2.2.1 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, Deff is assumed as a function of the particle-wall distance:

Deff=Deff(y)(14)

where y denotes the wall-distance in the normal direction in a boundary layer flow. Moreover, it is commonly accepted that Deff 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

Deff(y)=FD(y)Deff(15)

where Deff is a constant bulk diffusivity and FD(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

Deff(y)2cx2+y(Deff(y)cy+Deff(y)kTdϕdyc)u(y)cxλc+S=0(16)

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):

Deffcy=γc,y=0(17)

where

γ=Deff[0δϕ(FD(y)1eϕ(y)/kT1)dy]1(18)

FIGURE 1
www.frontiersin.org

FIGURE 1. Decoupling of the transport and deposition problems: solutions for the wall region (cw) and the bulk region (cw) 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.

The coefficient γ 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

Deffcn=γc(19)

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

cdt=λcd+γc(20)

where cd 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 γ, which has the physical dimensions of a velocity. Desorption mechanisms can be included as well by adding a corresponding term Eq 20,

cdt=(λ+δ)cd+γc(21)

with δ being the desorption rate constant (Wood et al., 2004).

2.3 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.

FIGURE 2
www.frontiersin.org

FIGURE 2. Parallel plates geometry: with respect to the flow, the x and y coordinates denote the longitudinal and transversal directions. Inlet and outlet boundaries are located at x=0 and x=L respectively, while wall boundaries are located at y=±H.

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.

2.3.1 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](22)

with

u(y)=um(1y2H2)(23)

where um is the maximum profile velocity and H is half the channel width as depicted in Figure 2.

2.3.2 Particle Concentration Equation

The boundary value problem then becomes, in explicit cartesian coordinates,

{D(2cx2+2cy2)u(y)cxλc+S=0(λ+δ)cd+γc=0,y=±HDcy=γcδcd,y=±H(24)

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:

x^=xL
y^=yH
u^(y^)=u(Hy^)um=1y^2
c^(x^,y^)=c(Lx^,Hy^)c0
c^d(x^)=cd(Lx^)Hc0
S^(x^,y^)=H2DS(Lx^,Hy^)c0

where c0 is a concentration value typical of the problem. The equations are then rewritten as

{f22c^x^2+2c^y^2fPeu^(y^)c^x^Dac^+S^=0c^d=ShDa+Δc^,y^=±1c^y^=DaShDa+Δc^,y^=±1(25)

where non-dimensional groups are defined as

Re=umHν(Reynolds number)
Sc=νD(Schmidt number)
Pe=ReSc=umHD(Péclet number)
Da=λH2D(Damköhler number)
Sh=γHD(Sherwood number)
Δ=δH2D(desorption number)
f=HL(aspect ratio)

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 Pef. 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

{2c^y^2fPeu^(y^)c^x^Dac^+S^=0c^d=ShDa+Δc^,y^=±1c^y^=DaShDa+Δc^,y^=1c^y^=0,y^=0c^=c^in=cinc0,x^=0(26)

where cin 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

S(x,y)=S¯S˜(x,y)(27)

where S¯ is a representative value of the source, e.g., its average value over the domain. The choice of c0 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:

c0=cin(28)

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

c0=S¯λ(29)

When decay is negligible, the reference concentration c0 can be chosen as

c0={H2S¯D,Sh1HS¯γ,Sh1(30)

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:

c^(x^,y^)=n=1Xn(x^)Yn(y^)(31)

2.3.3 Derivation of the Analytical Solution

The functions Yn(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:

n=1[fPeYndXndx^1u^(y^)(d2Yndy^2DaYn)Xn]=S^(x^,y^)u^(y^)(32)

The eigenvalue problem is therefore stated as

Yn=βn2Yn(33)

where

=1u^(y^)(d2dy^2Da)(34)

and βn2 being the eigenvalue associated to the eigenfunction Yn. The problem must be equipped with the corresponding boundary conditions:

{d2Yndy^2+(βn2Daβn2y^2)Yn=0dYndy^=DaShDa+ΔYn,y^=1dYndy^=0,y^=0(35)

It is easily verified that βn=0 is not an eigenvalue. Then, through the simple change of variable η=2βny^, Eq 35 can be recast into one form of the parabolic cylinder equation (DLMF, 2020, Ch. 12):

d2Yndη2[η24βn2(1Daβn2)]Yn=0(36)
Yn(y^)=A1,nY1,n(y^)+A2,nY2,n(y^)(37)

where

Y1,n(y^)=eβny^2/2M(1βn4+Da4βn,12;βny^2)(38)
Y2,n(y^)=2βny^eβny^2/2M(3βn4+Da4βn,32;βny^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:

M(a,b;z)=1+abz+a(a+1)b(b+1)z22!+=k=0(a)k(b)kzkk!(40)

where

(a)k=a(a+1)(a+k1)=s=1k(a+s1)(41)

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),

M(a,b;0)=1(42)
ddzM(a,b;z)=abM(a+1,b+1;z),(43)

it can be easily shown that the symmetry condition at y^=0 implies A2,n=0. The wall boundary condition requires

DaShDa+ΔY1,n(1)+Y1,n(1)=0(44)

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

Yn(y^)=Y1,n(y^)(45)

whose general solution (back in terms of y^) can be expressed as

The separated Eq 32 now reads

n=1(dXndx^+βn2fPeXn)Yn=1fPeS^(x^,y^)u^(y^)(46)

The eigenfunctions Yn are orthogonal with respect to the scalar product defined as

f,g=1+1u^(y^)fgdy^(47)

with

Ym,Yn=1+1u^(y^)YmYndy^=Cnδm,n(48)
Cn=Yn,Yn=1+1u^(y^)Yn2dy^(49)

This is exploited to transform Eq 46 into

dXndx^+βn2fPeXn=1fPeSn(x^)(50)

whose general solution is

Xn(x^)=exp(βn2fPex^)[Xn(0)+1fPe0x^exp(βn2fPeξ)Sn(ξ)dξ](51)

To obtain this last form, S^ has been expanded as

S^(x^,y^)=u^(y^)n=1Sn(x^)Yn(52)
Sn(x^)=1Cn1+1S^(x^,y^)Yndy^(53)

and the arbitrary constants Xn(0) have been determined from the remaining inlet boundary condition (expanded as well):

Xn(0)=1Cn1+1u^(y^)c^inYndy^=DnCnc^in(54)

where

Dn=1,Yn=1+1u^(y^)Yndy^(55)

The normalization constants Cn and Dn can be found by numerical integration of Eqs 49, 55, respectively. More details on their computation are reported in the Supplementary Appendix.

3 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

ShYn(1)+Yn(1)(56)

with

Yn(y^)=eβny^2/2M(1βn4+Da4βn,12;βny^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 36, 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.

FIGURE 3
www.frontiersin.org

FIGURE 3. Dependence of the 1st eigenvalue β1 on model parameters Sh and Da.

FIGURE 4
www.frontiersin.org

FIGURE 4. Dependence of the 2nd eigenvalue β5 on model parameters Sh and Da.

FIGURE 5
www.frontiersin.org

FIGURE 5. Dependence of the 5th eigenvalue β5 on model parameters Sh and Da.

FIGURE 6
www.frontiersin.org

FIGURE 6. Dependence of the 10th eigenvalue β10 on model parameters Sh and Da.

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.

TABLE 1
www.frontiersin.org

TABLE 1. Values of Da and Sh selected for verification cases.

TABLE 2
www.frontiersin.org

TABLE 2. Model parameters used in all verification cases.

To simplify the analysis, the inlet concentration is set to zero (cin=0). Furthermore, to allow for a direct comparison between different test cases, the reference concentration value has been selected as

c0=H2S¯D(58)

It follows that

S^(x^,y^)=S(Lx^,Hy^)S¯=S˜(Lx^,Hy^)(59)

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:

S˜(Lx^,Hy^)=π2cos(π2y^)(60)

Concentration profiles obtained for the nine test cases (Table 1) are shown in Figures 79.

FIGURE 7
www.frontiersin.org

FIGURE 7. Comparison of concentration profiles obtained with the proposed transport model implemented in OpenFOAM () and the corresponding analytical solutions () for cases 1, 2 and 3 (from top to bottom).

FIGURE 8
www.frontiersin.org

FIGURE 8. Comparison of concentration profiles obtained with the proposed transport model implemented in OpenFOAM () and the corresponding analytical solutions () for cases 4, 5 and 6 (from top to bottom).

FIGURE 9
www.frontiersin.org

FIGURE 9. Comparison of concentration profiles obtained with the proposed transport model implemented in OpenFOAM () and the corresponding analytical solutions () for cases 7, 8 and 9 (from top to bottom).

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 × 104 hexahedral volumes, with respectively 500 and 100 divisions on the longitudinal and transversal directions. All 9 cases shown similar convergence behavior, with approximately 5 × 104 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).

TABLE 3
www.frontiersin.org

TABLE 3. Computational times for the verification cases. 5 × 104 iterations are performed on 5 × 104 volumes with eight CPUs.

4 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.

Author Contributions

All authors contributed to conception and design of the study. AD performed the theoretical and technical analysis and wrote the first draft of the manuscript. SL, FG and AC supervised the analysis. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This project has received funding from the Euratom research and training programme 2014–2018 under grant agreement no. 847527.

Disclaimer

The content of this paper does not reflect the official opinion of the European Union. Responsibility for the information and/or views expressed therein lies entirely with the authors.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2021.692627/full#supplementary-material

References

Aufiero, M., Cammi, A., Geoffroy, O., Losa, M., Luzzi, L., Ricotti, M. E., et al. (2014). Development of an OpenFOAM Model for the Molten Salt Fast Reactor Transient Analysis. Chem. Eng. Sci. 111, 390–401. doi:10.1016/j.ces.2014.03.003

CrossRef Full Text | Google Scholar

Baes, C. F. (1974). The Chemistry and Thermodynamics of Molten Salt Reactor Fuels. J. Nucl. Mater. 51, 149–162. doi:10.1016/0022-3115(74)90124-X

CrossRef Full Text | Google Scholar

Balboa Usabiaga, F., Xie, X., Delgado-Buscalioni, R., and Donev, A. (2013). The Stokes-Einstein Relation at Moderate Schmidt Number. J. Chem. Phys. 139, 214113. doi:10.1063/1.4834696

PubMed Abstract | CrossRef Full Text | Google Scholar

Bowen, B. D., Levine, S., and Epstein, N. (1976). Fine Particle Deposition in Laminar Flow through Parallel-Plate and Cylindrical Channels. J. Colloid Interf. Sci. 54, 375–390. doi:10.1016/0021-9797(76)90317-9

CrossRef Full Text | Google Scholar

Brenner, H. (1961). The Slow Motion of a Sphere through a Viscous Fluid towards a Plane Surface. Chem. Eng. Sci. 16, 242–251. doi:10.1016/0009-2509(61)80035-3

CrossRef Full Text | Google Scholar

Cervi, E., Lorenzi, S., Cammi, A., and Luzzi, L. (2019a). Development of a Multiphysics Model for the Study of Fuel Compressibility Effects in the Molten Salt Fast Reactor. Chem. Eng. Sci. 193, 379–393. doi:10.1016/j.ces.2018.09.025

CrossRef Full Text | Google Scholar

Cervi, E., Lorenzi, S., Cammi, A., and Luzzi, L. (2019b). Development of an SP3 Neutron Transport Solver for the Analysis of the Molten Salt Fast Reactor. Nucl. Eng. Des. 346, 209–219. doi:10.1016/j.nucengdes.2019.03.001

CrossRef Full Text | Google Scholar

Compere, E., Kirslis, S., Bohlmann, E., Blankenship, F., and Grimes, W. (1975). “Fission Product Behavior in the Molten Salt Reactor experiment.”. Tech. Rep. ORNL–, 4865, 4077644 (Oak Ridge National Laboratory) Oak Ridge (TN). doi:10.2172/4077644

CrossRef Full Text | Google Scholar

Delpech, S., Merle-Lucotte, E., Heuer, D., Allibert, M., Ghetta, V., Le-Brun, C., et al. (2009). Reactor Physic and Reprocessing Scheme for Innovative Molten Salt Reactor System. J. Fluorine Chem. 130, 11–17. doi:10.1016/j.jfluchem.2008.07.009

CrossRef Full Text | Google Scholar

Fiorina, C., Hursin, M., and Pautz, A. (2017). Extension of the GeN-Foam Neutronic Solver to SP3 Analysis and Application to the CROCUS Experimental Reactor. Ann. Nucl. Energ. 101, 419–428. doi:10.1016/j.anucene.2016.11.042

CrossRef Full Text | Google Scholar

Fiorina, C., Kerkar, N., Mikityuk, K., Rubiolo, P., and Pautz, A. (2016). Development and Verification of the Neutron Diffusion Solver for the GeN-Foam Multi-Physics Platform. Ann. Nucl. Energ. 96, 212–222. doi:10.1016/j.anucene.2016.05.023

CrossRef Full Text | Google Scholar

DLMF (2020). “NIST Digital Library of Mathematical Functions. Release 1.1.0 of 2020-12-15,”. Editors F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clarket al.

Google Scholar

Grimes, W. R. (1970). Molten-salt Reactor Chemistry. Nucl. Appl. Tech. 8, 137–155. doi:10.13182/NT70-A28621

CrossRef Full Text | Google Scholar

Guha, A. (2008). Transport and Deposition of Particles in Turbulent and Laminar Flow. Annu. Rev. Fluid Mech. 40, 311–341. doi:10.1146/annurev.fluid.40.111406.102220

CrossRef Full Text | Google Scholar

Hébert, A. (2010). “Multigroup Neutron Transport and Diffusion Computations,” in Handbook of Nuclear Engineering. Editor D. G. Cacuci (Boston, MA: Springer US), 751–911. doi:10.1007/978-0-387-98149-9_8

CrossRef Full Text | Google Scholar

Issa, R. I. (1986). Solution of the Implicitly Discretised Fluid Flow Equations by Operator-Splitting. J. Comput. Phys. 62, 40–65. doi:10.1016/0021-9991(86)90099-9

CrossRef Full Text | Google Scholar

Johansson, F. (2017). Arb: Efficient Arbitrary-Precision Midpoint-Radius Interval Arithmetic. IEEE Trans. Comput. 66, 1281–1292. doi:10.1109/TC.2017.2690633

CrossRef Full Text | Google Scholar

Johansson, F. (2019). Computing Hypergeometric Functions Rigorously. ACM Trans. Math. Softw. 45, 1–26. doi:10.1145/3328732

CrossRef Full Text | Google Scholar

Kedl, R. J. (1972). “Migration of a Class of Fission Products (Noble Metals) in the Molten-Salt Reactor experiment,”. Tech. Rep. ORNL-TM–3884, 4471292(Oak Ridge National Laboratory) Oak Ridge (TN). doi:10.2172/4471292

CrossRef Full Text | Google Scholar

Launder, B. E., and Spalding, D. B. (1974). The Numerical Computation of Turbulent Flows. Comp. Methods Appl. Mech. Eng. 3, 269–289. doi:10.1016/0045-7825(74)90029-2

CrossRef Full Text | Google Scholar

Leppänen, J., Pusa, M., Viitanen, T., Valtavirta, V., and Kaltiaisenaho, T. (2015). The Serpent Monte Carlo Code: Status, Development and Applications in 2013. Ann. Nucl. Energ. 82, 142–150. doi:10.1016/j.anucene.2014.08.024

CrossRef Full Text | Google Scholar

Liu, B. Y. H., and Agarwal, J. K. (1974). Experimental Observation of Aerosol Deposition in Turbulent Flow. J. Aerosol Sci. 5, 145–155. doi:10.1016/0021-8502(74)90046-9

CrossRef Full Text | Google Scholar

Marino, A., Peltomäki, M., Lim, J., and Aerts, A. (2020). A Multi-Physics Computational Tool Based on CFD and GEM Chemical Equilibrium Solver for Modeling Coolant Chemistry in Nuclear Reactors. Prog. Nucl. Energ. 120, 103190. doi:10.1016/j.pnucene.2019.103190

CrossRef Full Text | Google Scholar

OpenFOAM (2021). The OpenFOAM Foundation. Online; (accessed 20 January. 2021).

Google Scholar

Patankar, S., and Spalding, D. (1972). A Calculation Procedure for Heat, Mass and Momentum Transfer in Three-Dimensional Parabolic Flows. Int. J. Heat Mass Transfer 15, 1787a–1806. doi:10.1016/0017-9310(72)90054-3

CrossRef Full Text | Google Scholar

Pearson, J. W., Olver, S., and Porter, M. A. (2017). Numerical Methods for the Computation of the Confluent and Gauss Hypergeometric Functions. Numer. Algorithms 74, 821–866. doi:10.1007/s11075-016-0173-0

CrossRef Full Text | Google Scholar

Price, T., Chvala, O., and Bereznai, G. (2020). A Dynamic Model of Xenon Behavior in the Molten Salt Reactor Experiment. Ann. Nucl. Energ. 144, 107535. doi:10.1016/j.anucene.2020.107535

CrossRef Full Text | Google Scholar

Prieve, D. C., and Ruckenstein, E. (1976). Rates of Deposition of Brownian Particles Calculated by Lumping Interaction Forces into a Boundary Condition. J. Colloid Interf. Sci. 57, 547–550. doi:10.1016/0021-9797(76)90232-0

CrossRef Full Text | Google Scholar

Ruckenstein, E., and Prieve, D. C. (1973). Rate of Deposition of Brownian Particles under the Action of london and Double-Layer Forces. J. Chem. Soc. Faraday Trans. 2 69, 1522–1536. doi:10.1039/F29736901522

CrossRef Full Text | Google Scholar

Serp, J., Allibert, M., Beneš, O., Delpech, S., Feynberg, O., Ghetta, V., et al. (2014). The Molten Salt Reactor (MSR) in Generation IV: Overview and Perspectives. Prog. Nucl. Energ. 77, 308–319. doi:10.1016/j.pnucene.2014.02.014

CrossRef Full Text | Google Scholar

Tiberga, M., de Oliveira, R. G. G., Cervi, E., Blanco, J. A., Lorenzi, S., Aufiero, M., et al. (2020). Results from a Multi-Physics Numerical Benchmark for Codes Dedicated to Molten Salt Fast Reactors. Ann. Nucl. Energ. 142, 107428. doi:10.1016/j.anucene.2020.107428

CrossRef Full Text | Google Scholar

Walker, S. A., and Ji, W. (2021). Species Transport Analysis of noble Metal Fission Product Transport, Deposition, and Extraction in the Molten Salt Reactor experiment. Ann. Nucl. Energ. 158, 108250. doi:10.1016/j.anucene.2021.108250

CrossRef Full Text | Google Scholar

Wood, B. D., Quintard, M., and Whitaker, S. (2004). Estimation of Adsorption Rate Coefficients Based on the Smoluchowski Equation. Chem. Eng. Sci. 59, 1905–1921. doi:10.1016/j.ces.2003.12.021

CrossRef Full Text | Google Scholar

Nomenclature

Latin symbols

Ar,g0 reference gth group, rth neutron reaction cross-section log. coeff.

c particle (vol.) concentration

c^ non-dim. particle (vol.) concentration

c^d non-dim. deposited particle (surf.) conc.

c^in non-dim. inlet particle (vol.) concentration

c0 ref. particle (vol.) concentration

cd deposited particle (surf.) concentration

ck kth family del. neutron prec. (vol.) concentration

cp specific heat capacity

cin inlet particle (vol.) concentration

D laminar particle diffusivity

dp particle diameter

Dt turbulent particle diffusivity

Deff effective particle diffusivity

Deff bulk effective particle diffusivity

Dn,g gth group neutron diffusion coefficient

Da Damköhler number

f channel aspect ratio

FD particle diffusivity correction factor

g gravitational acceleration

H channel half-width

kB Boltzmann constant

keff effective multiplication factor

L channel length

M confluent hypergeometric function

p pressure

Pe Péclet number

Pr Prandtl number

Prt turbulent Prandtl number

q vol. energy source

Re Reynolds number

S particle (vol.) source

S¯ average particle (vol.) source

S^ non-dim. particle (vol.) source

S˜ unit-average particle (vol.) source

Sc Schmidt number

Sct turbulent Schmidt number

Sh Sherwood number

T temperature

T0 reference temperature (buoyancy)

T0Σ reference temperature (cross-sections)

u longitudinal velocity comp. velocity

u^ non-dim. longitudinal velocity comp.

u longitudinal velocity comp. velocity

um maximum profile velocity

vg avg. gth group neutron velocity

x longitudinal coordinate

x^ non-dim. longitudinal coord.

y transversal coord.

y^ non-dim. transversal coord.

yf fission yield of transported fission product

Greek symbols

α laminar thermal diffusivity

αt turbulent thermal diffusivity

αeff effective thermal diffusivity

βd total delayed neutrons fraction (kβd,k)

βd,k kth family delayed neutrons fraction

βn (square root of) nth eigenvalue

βT vol. thermal expansion coeff.

Δ desorption number

δ desorption rate constant

δϕ wall region thickness

γ deposition velocity

λ particle decay constant

λk kth family del. neutron prec. decay constant

ρ fluid density

Σa,g gth group absorption cross-section

Σf,g gth group fission cross-section

Σr,g gth group, rth neutron reaction cross-section

Σr,g0 reference gth group, rth neutron reaction cross-section

Σs,gh g-to-h-th group scattering cross-section

χd,g gth energy group delayed neutron spectrum

χp,g gth energy group prompt neutron spectrum

ν¯g avg. gth group neutrons emitted per fission

ν laminar kinematic viscosity

νeff effective kinematic viscosity

νt turbulent kinematic viscosity

ϕ particle-wall interaction potential

φg gth group integrated neutron flux

Keywords: molten salt fast reactor, molten salt, fission products, particle deposition, multiphysics

Citation: Di Ronco A, Lorenzi S, Giacobbo F and Cammi A (2021) An Eulerian Single-Phase Transport Model for Solid Fission Products in the Molten Salt Fast Reactor: Development of an Analytical Solution for Verification Purposes. Front. Energy Res. 9:692627. doi: 10.3389/fenrg.2021.692627

Received: 08 April 2021; Accepted: 01 June 2021;
Published: 29 June 2021.

Edited by:

Wenzhong Zhou, Sun Yat-Sen University, China

Reviewed by:

Yu Ma, Sun Yat-Sen University, China
Jiankai Yu, Massachusetts Institute of Technology, United States

Copyright © 2021 Di Ronco, Lorenzi, Giacobbo and Cammi. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Antonio Cammi, antonio.cammi@polimi.it

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.