Development and test of a novel verification scheme applied to the neutronic modelling of Molten Salt Reactors

This paper presents the extension of a method to verify transient neutron transport solvers earlier developed for reactors with non-moving fuel, to the case of Molten Salt Reactors (MSRs). This method is based on the extraction of the point-kinetic response of a nuclear reactor excited by a mono-chromatic perturbation and on its subsequent comparison with its expected functional dependence. Whereas a simple expression for this dependence exists for systems with ﬁxed fuel, this is not the case for MSRs, as high-lighted in many past studies. A workaround is nevertheless proposed in this work, thus giving the possibility to use a similar veriﬁcation method to the case of MSRs. The method is applied to a simple dynamic MSR solver, demonstrating the capabilities of the technique. Contrary to other veriﬁcation methods for which the system has to be simpliﬁed so that analytical solutions can be derived, the present method can be applied to any heterogeneous system.


Introduction
Molten Salt Reactors (MSRs) represent, from a reactor physics and fuel processing viewpoint in particular, systems that significantly differ from today's reactors.The difference arises from the movement and possible online reprocessing of the fuel in MSRs, as opposed to static fuel in present reactors.Because of the interplay between the prompt and delayed neutrons occurring at different time scales, the movement of the fuel gives rise to peculiarities in MSRs, both in steady-state and non-steady-state conditions.Such specificities have been known since the early developments of MSRs (Prince et al., 1968;Haubenreich, 1969;Haubenreich and Engel, 1970).
In the early 2000s, the MSR concept was selected by the Generation IV International Forum as one of the six Generation IV reactor technologies to be further developed.The renewed interest in MSRs resulted in many research groups developing modelling capabilities specifically targeted at MSRs.As part of the software development, efforts are dedicated to the verification and validation of the tools -see e.g., (Tiberga et al., 2020).
Typically, verification relies on the demonstration that the numerical implementation of a chosen model is correct.Validation, on the other hand, is based on the successful comparison between the results of the solver being considered and either experimental data or the results of other already verified/validated solvers (Boehm, 1981;Oberkampf and Trucano, 2008;Oberkampf et al., 2004).Although such solvers might use other approximations, their solutions can be considered as the reference solution.In the area of code verification, a simple enough system configuration is chosen, so that an analytical solution can be calculated, to which the numerical solver can be compared, while refining the discretization of the variables.Alternatively, another solver that has already been verified and that relies on the same approximations of the representation of the physical phenomena being modelling can be used.In recent years, the Method of Manufactured Solutions (MMS) has also received a particular attention for code verification (Roache et al., 2019).In this method, a mathematical expression is chosen.Putting this expression in the non-discretized balance equations allows computing a source term to those.This source term can then be inserted into the discretized equations in the code being verified.When refining the discretization of the variables, the code should thus return the artificially manufactured original mathematical expression.
This paper presents the extension of a validation/verification method of transient neutron transport solvers earlier developed in (Demazière et al., 2017) for the case of reactors with fixed fuel, to the case of MSRs.As detailed in (Demazière et al., 2017), the technique falls in the category of both code verification and validation.This method is based on the extraction of the point-kinetic response of a system excited by a monochromatic perturbation as computed by a modelling software and on its subsequent comparison with its expected functional dependence when the frequency of the perturbation is varied.This functional dependence reduces to a simple analytical form for systems with fixed fuel.The point-kinetic response is then given by the well-known zero-power reactor transfer function multiplied by the reactivity effect of the applied perturbation (Bell and Glasstone, 1970).In the case of MSRs, the formulation of an equivalent zero-power reactor transfer function is much more involved, and no simple expression exists between the reactivity effect and the corresponding point-kinetic response, as pointed out in several past studies -see e.g., (Pázsit et al., 2014).The main reason behind the impossibility to write a closed form of the zero-power reactor transfer function in the case of MSRs lies with the spatial dependence of the adjoint functions being different with respect to the neutron flux and to the precursors of delayed neutrons, respectively.In the case of a reactor with fixed fuel, these two adjoint functions have the same spatial dependence.In the case of MSRs, the spatial dependence of the adjoint function associated to the precursors of delayed neutrons is fundamentally different from the one associated to the neutron flux, due to the transport of the precursors with the moving fuel.
Despite this difference and the complications arising from it, the work reported hereafter presents a workaround.The novel method relies on the extraction of the point-kinetic component of both the neutron flux and the precursors of delayed neutrons, respectively, and on the estimation of corresponding variations of the so-called shape functions, for both the neutron flux and the precursors of delayed neutrons.As a function of the frequency of the applied perturbation, functional dependences between those four quantities exist and need to be verified.The inability of a code output to reproduce such functional dependencies might correspond to a deficiency of the software.The method reported in the following thus offers the possibility to verify the proper implementation of a neutron transport solver for the case of MSRs without turning to the MMS or analytical/semi-analytical solutions, for which simplified homogeneous or piece-wise homogeneous systems in one or two dimensions at the most need to be considered.The novel method can instead be applied to three-dimensional heterogeneous systems.
This paper is structured as follows.The next section defines the system being considered, as well as the modelling framework used for deriving the verification framework.As explained in this paper, the verification method relies on the excitation of the system by a mono-chromatic perturbation expressed as fluctuations of the macroscopic cross-sections, as compared to a stationary state of the reactor.Both the stationary reactor conditions and the deviations from the corresponding steady-state conditions need to be considered.The former is presented in the third section of this paper, together with the formulation of the adjoint functions necessary for the verification method.The latter is considered in the fourth section of the paper, assuming a factorization of the neutron flux and the precursors of delayed neutrons into an amplitude factor and a shape function.Such factorizations are required for deriving the linearized point-kinetic equations, presented in the fifth section of the paper.The verification scheme is then introduced in the sixth section and the demonstration of the method in the seventh section.The paper ends with some conclusions and remarks on the applicability of the method.Although the method is derived in the most general case of three-dimensional systems and could be readily used for verifying any three-dimensional solver, the method is applied to and demonstrated on a one-dimensional system for the sake of simplicity and illustration.

Definition of the system and of the associated modelling framework
As indicated in the introduction, the verification method is based on the estimation of the so-called linearized neutron noise in the frequency domain, which will be formally defined in the subsequent sections of this paper, and on the extraction of the point-kinetic components for the neutron noise and the noise in the concentration of the precursors of delayed neutrons.Pioneering work related to neutron noise in the specific case of MSRs is reported in depth in (Pázsit and Dykin, 2017).
Although the verification method reported hereafter can be applied to neutron transport solvers, it will be demonstrated in the framework of two energy group diffusion theory, since most of the neutronic core simulators are based on such a resolution in energy and angle.The methodology followed to derive the governing equations describing the noise in the neutron density and in the concentration of precursors of delayed neutrons is standard and largely inspired from past studies in this area (Pázsit and Dykin, 2017;Pázsit and Jonsson, 2011;Jonsson and Pázsit, 2011;Pázsit et al., 2012;Dykin and Pázsit, 2016).In this work, threedimensional heterogeneous systems are however considered, whereas one-dimensional homogeneous systems were investigated in (Pázsit and Dykin, 2017;Pázsit and Jonsson, 2011;Jonsson and Pázsit, 2011;Pázsit et al., 2012;Dykin and Pázsit, 2016).Moreover, the velocity of the fuel is space-dependent, as opposed to a space-independent velocity assumed in the models presented in (Pázsit and Dykin, 2017;Pázsit and Jonsson, 2011;Jonsson and Pázsit, 2011;Pázsit et al., 2012;Dykin and Pázsit, 2016).The non-homogeneous nature of the velocity has some implication for the estimation of the adjoint operators, to which we will return in the next Section.Finally, the last major difference between the model used in this paper and earlier investigations lies with the boundary condition used in such investigations.Instead of a simple decay of precursors outside of the core, a more generic inlet boundary condition in the concentration of the neutron precursors is implemented.Such a choice is more suitable for modelling a three-dimensional system with on-line reprocessing of the fuel.This will have some implication for the estimation of the concentration of the precursors of delayed neutrons and on their corresponding adjoint functions.
For a three-dimensional heterogeneous MSR of cylindrical shape in which the fuel velocity is assumed to be timeindependent and mono-directional along the vertical upward axis, i.e., with z being the unit upward vector, the time-and spacedependent balance equations in two group diffusion theory and one group of delayed neutrons read as: In the above equation: / g r; t ð Þ represents the space-and time-dependent scalar neutron flux in the energy group g.C r; t ð Þ represents the space-and time-dependent concentration of the precursors of delayed neutrons.v g;0 r ð Þ is the neutron speed, assumed to be only spacedependent.k 0 r ð Þ is the decay constant of the precursors of delayed neutrons, assumed to be only space-dependent.b 0 r ð Þ is the fraction of delayed neutrons, assumed to be only space-dependent.R 1;1 r; t ð Þ is a space-and time-dependent operator acting upon / 1 r; t ð Þ as: is the diffusion coefficient in the energy group g, assumed to be only space-dependent.tR f ;g r; t ð Þ is the space-and time-dependent product between the average number of neutrons released by fission event and the macroscopic fission cross-section, in the energy group g.R a;g r; t ð Þ is the space-and time-dependent macroscopic absorption cross-section in the energy group g.R r r; t ð Þ is the space-and time-dependent macroscopic removal cross-section.I r; t ð Þ is a time-and space-dependent operator acting upon C r; t ð Þ as: which, in the present case, further simplifies into: The governing equations given by Eq. ( 2) are complemented by the following boundary conditions: where r b represents the boundary of the system, r rad the radial only boundary, and r in the inlet only boundary.r rad and r in are thus subsets of r b .All distances related to the boundaries are assumed to be extrapolated.
As compared to earlier work, e.g., (Pázsit and Dykin, 2017;Pázsit and Jonsson, 2011;Jonsson and Pázsit, 2011;Pázsit et al., 2012;Dykin and Pázsit, 2016), the inlet boundary condition with respect to the precursors of delayed neutrons defined in Eq. ( 9) offers a larger flexibility for representing the possible online reprocessing of the moving fuel, the inherent mixing of the fuel between the outlet and the inlet, the possible removal and addition of precursors, etc. Specifying a given relationship between the outlet and inlet of the core with respect to the precursors of delayed neutrons, for instance a pure decay of the form C r in ; t (Pázsit and Dykin, 2017;Pázsit and Jonsson, 2011;Jonsson and Pázsit, 2011;Pázsit et al., 2012;Dykin and Pázsit, 2016) (with s representing the transit time in the external loop), can thus be considered as a subcase of the model considered in this work.

Definition of the static problem
The static problem is formulated by removing the timedependence in Eq. ( 2).This results in the following balance equations: where the subscript 0 represents time-independent quantities.When deriving the point-kinetic equations in Section 5, the adjoint functions associated with the neutron flux and the precursors of delayed neutrons will be needed.We introduce hereafter adjoint operators and functions, following common practice in reactor dynamics (Bell and Glasstone, 1970).For any spacedependent operator @ r ð Þ acting upon a space-dependent function f r ð Þ being a column vector having for components f 1 r ð Þ and f 2 r ð Þ for the fast energy and thermal energy groups, respectively, we thus define the corresponding adjoint operator @ þ r ð Þ and adjoint func- for the fast energy and thermal energy groups, respectively.The adjoint operator and functions should fulfil the following condition: with the inner product between two space-dependent functions u r ð Þ and w r ð Þ being defined, in a two-energy group structure, as: In the equation above, T denotes the transpose operator, and V represents the volume of the reactor core.Finding the adjoint operator corresponding to operators involving spatial derivatives requires some extra attention.In the balance equations given by Eq. ( 10), two types of operators involving spatial derivatives , one could demonstrate, using a double integration by parts, Gauss divergence theorem and the boundary conditions for the static neutron fluxes defined by Eq. ( 9), that: if the following boundary condition for the adjoint function of the neutron flux is introduced: For the operator $ Á u 0 r , one could demonstrate, using an integration by parts, Gauss divergence theorem and the boundary conditions for the precursors of delayed neutrons defined by Eq. ( 9), that: if the following boundary condition for the adjoint function of the neutron flux is introduced: and using cylindrical coordinates to describe the position in the reactor core.
It could further be demonstrated that the property of the adjoint defined in Eq. ( 11) requires the adjoint functions to fulfil the following balance equations: where Eqs. ( 13) and ( 15) were used.In Eq. ( 18), one has: which, in the present case, simplifies into: 4. Definition of the dynamic problem in its factorized form In order to later derive the point-kinetic approximation, the time-and space-dependent scalar neutron fluxes and concentration of the precursors of delayed neutrons are factorized as: with P t ð Þ and Q t ð Þ being the time-dependent amplitude factors associated with the scalar neutron fluxes and concentration of the precursors, respectively, whereas w g r; t ð Þ; g ¼ 1; 2 and u r; t ð Þ are the time-and space-dependent shape functions associated with the scalar neutron fluxes and concentration of the precursors, respectively.
The system is assumed to be at steady-state conditions until the time t ¼ 0, which represents the time at which a perturbation to the system is applied.For the sake of simplicity, we normalize the amplitude factors at time t ¼ 0 to unity as: This results in: and The boundary conditions expressed in Eq. ( 9) lead to: In order to be able to obtain the point-kinetic approximation (later derived in this Section), the following normalization condi-tions of the shape functions are introduced for the neutron flux (Bell and Glasstone, 1970): and for the delayed neutron precursors (Bell and Glasstone, 1970): where Eqs. ( 24) and ( 25) were used for the last equalities in Eqs. ( 28) and (30), respectively.Eqs. ( 28) and ( 30) also express the fact that the quantities time-independent.It should be mentioned that the conditions expressed by Eqs. ( 27) and ( 29) are the only possible normalizations if one wants to derive the point-kinetic equations, according to which balance equations for the time-dependent amplitude factors only are obtained.Any other normalization would result in time-derivatives associated with the shape functions remaining as well.
The three time-dependent equations given by Eq. ( 2) are then multiplied by , respectively, and integrated on the entire reactor volume.The adjoint balance equations given by Eq. ( 18) are multiplied by / 1 r; t ð Þ, / 2 r; t ð Þ, and C r; t ð Þ, respectively, taken in their factorized form (i.e., expressed using Eqs.( 21) and ( 22)), and integrated on the entire reactor volume.Taking the difference between these two sets of equations results in: and In the equations above, the various quantities are defined as: In the equations above, the fluctuations of any time-dependent and most generally space-dependent quantity X r; t ð Þ are defined as: with X 0 r ð Þ representing the steady-state value of X r; t ð Þ. Eqs. ( 31) and (32) were obtained making use of the fact that the steady-state values of all quantities appearing in the balance equations fulfil Eq. ( 10).Moreover, the following relationships were utilized: Eqs. ( 41) and ( 42) can be derived in the same manner as Eqs.( 13) and ( 15), using the boundary conditions of the shape functions introduced in Eq. ( 26), and more precisely w g r b ; t ð Þ¼0; g ¼ 1; 2; 8t and u r rad ; t ð Þ¼0; 8t.
Finally, the normalization conditions corresponding to Eqs. ( 27)-(30) were introduced to eliminate the volume integrals involving partial time-derivatives related to any of the shape functions.The normalization conditions are precisely introduced to eliminate such terms.It is only thanks to the normalization conditions defined by Eqs. ( 27) and ( 29) that expressions for the time variations of the amplitude factors P t ð Þ and Q t ð Þ alone are obtained, as earlier mentioned.Eqs. ( 31) and ( 32) correspond to the classical point-kinetic equations, although their expressions are more involved in the present case of MSRs as compared to the case of reactors with fixed fuel.These equations are comparable to the ones obtained elsewhere (see e.g.(Pázsit et al., 2014) in one-group diffusion theory and (Dykin and Pázsit, 2016) in two-group diffusion theory), if the same modelling assumptions are used.
Comparing the point-kinetics equations for MSRs with the ones for reactors with fixed fuel, beyond the fact that the shape functions are different with respect to the neutron fluxes and precursors of delayed neutrons, respectively, a new term, S t ð Þ, exists in the time-dependent equation associated to the amplitude factor of the precursors of delayed neutrons, as earlier pointed out by Pázsit et al. in (Pázsit et al., 2014).As Eq. ( 39) demonstrates, this term is not equal to zero in an MSR.This is in clear contrast with systems with fixed fuel, for which strict equality with zero would be obtained, since in such systems the concentration of the neutron precursors has the same boundary conditions as for the scalar neutron flux.The fact that such a simple boundary condition is not applicable in MSRs makes such systems more challenging when deriving reactor dynamics approximations.

Linearization of the point-kinetic equations
The normalization of the shape functions given by Eqs. ( 28) and (30), respectively, leads to: Then Eqs. ( 31) and ( 32) can be rewritten as: Following the same approach as in (Pázsit et al., 2014), all timedependent quantities are expressed as sums between their steadystate values and their deviations from steady-state -see Eq. ( 40).Neglecting second-order perturbation terms and taking a Fourier transform gives: which represent the linearized point-kinetic equations in the frequency domain.For the sake of simplicity, the same notations were used for functions f t ð Þ depending on time t and their Fouriertransform counterpart f x ð Þ, with x being the angular frequency.
In first order, the various terms appearing in Eq. ( 46) are defined as: where the normalization conditions given by Eqs. ( 28) and (30) were used, as well as Eq. ( 17).A closer examination of Eqs. ( 48), ( 49) and ( 51) reveals that the knowledge of the fluctuations of the shape functions dw g r; The fluctuations of the shape functions can be evaluated as detailed below.
Using a first-order approximation, one has for the neutron flux: The normalization condition given by Eq. ( 27), together with the initial condition dw g r; t ¼ 0 ð Þ¼0; g ¼ 1; 2 derived from Eq. ( 24), results in: which gives in the frequency domain: From Eqs. ( 54) and ( 56), one concludes that: Once dP x ð Þ has been evaluated from Eq. ( 58), the fluctuations of the shape functions are thus given, via Eq.( 54), as: Likewise, using a first-order approximation, one has for the precursors of delayed neutrons: The normalization condition given by Eq. ( 29), together with the initial condition du r; t ¼ 0 ð Þ¼0 derived from Eq. ( 25), results in: which gives in the frequency domain: From Eq. ( 62), one concludes that: Once dQ x ð Þ has been evaluated from Eq. ( 64), the fluctuations of the shape functions are thus given, using Eq. ( 60), as: 6. Derivation of the verification method Eq. ( 46) can be compared to the one obtained in the case of reactors with fixed fuel, for which the following relation holds (Bell and Glasstone, 1970): with the zero-power reactor transfer function G 0 x ð Þ being defined as: In the case of reactors with fixed fuel, a simple analytical expression exists between the point-kinetic component of the induced neutron noise dP x ð Þj u z;0 ¼0 and the reactivity noise dq x ð Þ, via the zero-power reactor transfer function G 0 x ð Þ.Once the perturbation of the system is defined, dq x ð Þ can be estimated accordingly.Multiplying dq x ð Þ by G 0 x ð Þ thus gives the point-kinetic component of the induced neutron noise, which could also be compared to the point-kinetic component extracted from the full solution to the problem via Eq.( 58).Varying the angular frequency of the applied perturbation allows verifying whether the extracted point-kinetic response from a code output coincides with its expected analytical solution.This is exactly the principle of the verification method proposed in (Demazière et al., 2017) for the case of reactors with fixed fuel.
In the case of MSRs, no simple analytical expression for the zero-power reactor transfer function exists.Nevertheless, a func- and dS x ð Þ exists, as Eq. ( 46) highlights.Since the last four terms themselves only depend on the deviation from point-kinetics of the neutron flux, dw g r; A verification method equivalent to the case of reactors with fixed fuel can thus be derived.The method is aimed at demonstrating that the functional dependencies f and g are verified between Assuming that the neutron noise d/ g r; x ð Þ; g ¼ 1; 2 and the noise dC r; x ð Þ in the concentration of the precursors of delayed neutrons induced by a given perturbation are estimated by a software, the verification method proposed in this paper goes as follows: 1. dP x ð Þ and dQ x ð Þ are extracted from the code output using Eqs.( 58) and (65), respectively.

dP x
ð Þ and dQ x ð Þ are also estimated from Eq. ( 46) (with dP x ð Þ evaluated from the first equation and used on the right-hand side of the second equation for dQ x ð Þ).

dP x
ð Þ and dQ x ð Þ estimated at step 1 are compared with those estimated at step 4.
The process above is repeated for various angular frequencies x.
Step 5 then allows determining whether the software works as intended.

Application and demonstration of the verification method
For the purpose of illustration, the method derived above is applied to a simple MSR solver.The system being considered corresponds to the test case used in (Pázsit et al., 2014): a homogeneous one-dimensional reactor of size 300 cm, in which the fuel is circulating through the core at a homogeneous speed u z;0 of 50 cm/s.One energy group and one group of delayed neutrons are assumed.The reactor composition is thus entirely defined by its macroscopic data tR f ;0 ¼ 1:0075 Â 10 À2 cm À1 , R a;0 ¼ 10 À2 cm À1 , and D 0 ¼ 0:33 cm, and by its kinetic data b 0 ¼ 650 pcm and k 0 ¼ 0:1 s À1 .The fission cross-section is also adjusted in order to obtain a critical system.The dynamic problem is defined by introducing a mono-chromatic point-wise perturbation in the macroscopic absorption cross-section at the center of the reactor.A frequency-independent noise source is assumed (white noise).A time delay s of 8 s between the outlet and the inlet is assumed for the decay of the precursors of delayed neutrons, as: A one-energy group diffusion-based solver was developed to estimate the static solution and the fluctuations in neutron flux and in the concentration of the precursors of delayed neutrons in the frequency-domain directly.The adjoint functions to the static problem were then determined from the static forward solution following the derivation in (Pázsit et al., 2014), which yields: The spatial discretization of the solver is based on finite differences.In the present test case, the mesh size was set to 1 cm, with the discretized quantities evaluated at the faces of each mesh (point-scheme).
The verification process described in Section 6 is illustrated for the solver under investigation.Fig. 1 and Fig. 2 show the amplitude and phase, respectively, of the point-kinetic component dP x ð Þ of the neutron flux, and Fig. 3 and Fig. 4 show the amplitude and phase, respectively, of the point-kinetic component dQ x ð Þ of the precursors of delayed neutrons.The figures also include the relative differences computed, for the amplitude factor dP x ð Þ related to the neutron flux, as: and for the amplitude factor dQ x ð Þ related to the precursors of delayed neutrons, as: amplitude of dQ x ð Þ À amplitude of g dw 1 r; The differences in absolute terms are computed, for the phase of the amplitude factor dP x ð Þ related to the neutron flux, as: and for the phase of the amplitude factor dQ x ð Þ related to the precursors of delayed neutrons, as: phase of dQ x ð Þ À phase of g dw 1 r; In Eqs. ( 71)-( 74), dP x ð Þ and dQ x ð Þ are directly computed from the code output using Eqs.( 58) and ( 64), respectively, whereas the functions f and g, respectively, represent the evaluation of the corresponding terms dP x ð Þ and dQ x ð Þ using Eq. ( 46).As can be seen in those figures, the agreement between the point-kinetic components directly computed from Eqs. ( 58) and (64) for the neutron flux and the precursors of delayed neutrons, respectively, and computed using the expressions given by Eq. ( 46) is excellent, both regarding the amplitude and phase.Despite varying over several orders of magnitude as a function of frequency, the amplitudes of the point-kinetic terms are correctly predicted.The peaks present at given frequencies in the amplitude of dP x ð Þ and dQ x ð Þ correspond to resonances introduced by the recirculation of the precursors of delayed neutrons, as explained in, e.g., (Pázsit et al., 2014).The agreement between dP x ð Þ and dQ x ð Þ directly computed from the code output and the functions f and g, respectively, demonstrates that the solver is working satisfactorily.In case of discrepancies, the examination of those may help to identify possible ''bugs" in the solver.

Conclusions
In this paper, the verification method earlier developed in (Demazière et al., 2017) for the case of reactors with fixed fuel was extended to the case of reactors with moving fuel.Whereas in the former case, the extraction of the point-kinetic component from the results of a code can be directly compared with an analytical expression, the comparison in the latter case requires additional quantities, which can only be estimated from the simulations.Such quantities appear on the right hand-side of Eq. ( 46) for dP x ð Þ as: If the above was negligible in comparison with:

Fig. 4 .
Fig. 4. Variation of the phase of the point-kinetic component dQ x ð Þ of the precursors of delayed neutrons.

Fig. 1 .
Fig. 1.Variation of the amplitude of the point-kinetic component dP x ð Þ of the neutron flux.

Fig. 2 .Fig. 3 .
Fig. 2. Variation of the phase of the point-kinetic component dP x ð Þ of the neutron flux.