APPLICATION OF RESPONSE MATRIX METHOD TO TRANSIENT SIMULATIONS OF NUCLEAR SYSTEMS

Until recently, reactor transient problems were exclusively solved by approximate deterministic methods. The increase in available computing power made it feasible to approach the transient analyses with time-dependent Monte Carlo methods. These methods offer the first-principle solution to the space-time evolution of reactor power by explicitly tracking prompt neutrons, precursors of delayed neutrons and delayed neutrons in time and space. Nevertheless, a very significant computing cost is associated with such methods. The general benefits of the Monte Carlo approach may be retained at a reduced computing cost by applying a hybrid stochastic-deterministic computing scheme. Among such schemes are those based on the fission matrix and the response matrix formalisms. These schemes aim at estimating a variant of the Greens function during a Monte Carlo transport calculation, which is later used to formulate a deterministic approach to solving a space-time dependent problem. In this contribution, we provide an overview of the time-dependent response matrix method, which describes a system by a set of response functions. We have recently suggested an approach where the functions are determined during a Monte Carlo criticality calculation and are then used to deterministically solve the space-time behaviour of the system. Here, we compare the time-dependent response matrix solution with the transient fission matrix and the time-dependent Monte Carlo solutions for a control rod movement problem in a mini-core reactor geometry. The response matrix formalism results in a set of loosely connected equations which offers favourable scaling properties compared to the methods based on the fission matrix formalism.


INTRODUCTION
High computing cost of time-dependent Monte Carlo calculations is considered as the main limitation for a wide-spread application of such calculations. Computing times ranging from a couple of hundreds of CPU hours (for small kinetics problems in simple geometries) to multiple thousands of CPU hours (for problems in assembly and mini-core geometries) were reported [1,2]. Such computing times would extrapolate to multiple CPU years for realistic large-scale problems with feedbacks. Using a hybrid stochastic-deterministic method may allow retaining the benefits of a Monte Carlo calculation, such as accuracy and flexibility while reducing the computing cost. Alternatively, a space-time dependent Monte Carlo based solver could potentially be applied as a tool in time-dependent Monte Carlo calculations for improving the figure-of-merit of such calculations.
One such approach, called the Transient Fission Matrix (TFM) method was recently proposed [3]. The method describes the space-time behaviour of the system by a set of fission and time matrices which can be calculated by Monte Carlo codes. The matrices are then used in formulating kinetics equations for neutron and precursor populations, which can then be solved deterministically. The TFM method was demonstrated on a variety of reactor physics problems [3,4].
Recently, we proposed a similar approach based on the Response Matrix formalism [5]. Similarly to the fission matrix methods, the RMM also superimposes a spatial mesh over the system geometry and calculates neutron propagation rates using the mesh. However, unlike the fission matrix methods, the RMM computes a number of parameters (response functions) for each node. The response functions computed for a given node represent the probabilities with which neutrons entering a node through a surface, or born in the node, propagate through the node to other surfaces of the same node. This way, each node is coupled only to the neighbouring nodes, differently from the approaches using the fission matrix, where each node is coupled to all other nodes in the system.
In this contribution, we compare the three methods for a control rod movement problem in a minicore reactor geometry. We briefly introduce the methods in Section 2, then describe the test case in Section 3. The results are presented in Section 4. Finally we conclude in Section 5.

Time-Dependent Monte Carlo Method
In Time-Dependent Monte Carlo (TDMC) methods, the time variable is treated explicitly by simply assigning a time parameter (the "internal clock") for each transported particle. The clock is set to zero at the start of the simulation and is progressively updated while simulating the particle's history. If the particle crosses the time interval boundary, it is stopped exactly at the boundary and a new path length (and the corresponding flight time) is sampled with (possibly) updated crosssections. The total calculation time is split into time intervals, which are used to score the quantities of interest, introduce system feedbacks, perform population control, or adjust other calculation parameters (e.g. weight windows) [1]. The time intervals can be non-uniform, and different interval sizes can be used for different tasks, e.g. quantity scoring or population control [2]. Differently from deterministic methods, where the size of the time step can numerically influence the accuracy of the simulation results, the sizes of the time intervals in a TDMC calculation only determine the sizes of the tally bins and therefore the resolution in which the tallied quantity is represented [1].
Simulating nuclear reactor transients by a Monte Carlo process poses several challenges due to delayed neutrons and fission chain branching [1]. The large time-scale differences between prompt neutron chains and delayed neutron emissions yield considerable time spread between the subsequent delayed neutron-induced fission chains, which in turn results in high variance of the results. Branching of neutron chains during fission causes a significant spread in the neutron chain lengths and consequently yields high variance of the results as well. High variance in TDMC simulations means that a large number of samples has to be collected for acceptable statistical accuracy. This results in prohibitive computing time.

Transient Fission Matrix Method
The Transient Fission Matrix method relies on four fission matrices, weighted on spectra and multiplicities of prompt and delayed neutrons G χp,νp , G χp,ν d , G χ d ,νp , G χ d ,ν d , and a time matrix T χp,νp which stores the fission-to-fission time. Using the matrices, balance equations are written for the neutron N and precursor P f vectors as [3,4] where β f is the normalized delayed neutron fraction ( β f = 1), and λ f is the decay constant for precursor group f . Here, a single, effective fission-to-fission time l ef f is used to characterize the prompt time response. The l ef f is calculated using the time matrix as where (G χp,νp · T χp,νp ) is the matrix element-to-element multiplication, which gives the neutron production associated to the response time from node j to i, N p is the eigenvector of G χp,νp , and N * p is the eigenvector of a transposed fission matrix G T χp,νp .
Using the TFM method, the four fission matrices and the time matrix are scored during Monte Carlo criticality calculations, corresponding to different states of the system, or by using a correlated sampling approach during a single criticality calculation [4]. The time-dependent solution is then obtained by integrating Eqs. (1)-(2).

Time-Dependent Response Matrix Method
The response matrix method is based on the decomposition of a global spatial domain V into a set of N disjoint subvolumes (nodes) V i that satisfy V = V 1 ∪ V 2 ∪ · · · ∪ V N . Local particle balance is represented by relating partial currents incident on and generated in a node to the partial currents leaving from or absorbed in that node via a set of coefficients usually called response functions. Global particle balance is then defined by a set of equations relating the outgoing partial currents from one node to the partial currents incident on another node. Denoting the source vector as s p (t), the partial (incident) current vector as j − (t), and the precursor vector as c(t), the time-dependent response equations for these quantities can be formally expressed as [5] s where M is a "connectivity" or "topological" matrix that directs the outgoing current from one node as the incident current to another. Here, various R(t → t) terms stand for the response functions, which represent the expected contribution to the prompt source, delayed source, or partial currents at time t from an average prompt or delayed neutron born in a node, or a unitary current incident on the node surface, at time t . The subscripts indicate the origin of neutrons and include p for the prompt neutron source, d for the delayed neutron source, and c for the partial current. The superscript indicates the class of neutron produced. When deriving Eqs. (4)-(6) we assumed that the response functions can be evaluated during a Monte Carlo criticality calculation. Such evaluation of the Green's function drives the Transient Fission Matrix method [3,4] and the use of the response matrix method for acceleration of fission-source convergence in Monte Carlo criticality calculations [6].
Several models based on orthogonal expansion were previously proposed to treat the time integrals in Eqs. (4)-(6) [7,8]. Here, we choose to apply a simple approach originally proposed in the work of Sicilian and Prior [9]. First, a set of discrete times t 0 = 0, t 1 , t 2 , . . . , t n , at which the state variables s p (t), j − (t) and c(t) will be evaluated, is selected, with s p,n ≡ s p (t n ), j − n ≡ j − (t n ) and c n ≡ c(t n ). Next, the partial currents, j − (t), and the source rates s p (t) are assumed to vary linearly between the time points, while the precursor concentrations c(t) are assumed to remain constant between the successive time points. Applying the outlined assumptions on the time variable, we obtain the discretized equations c n+1,j = e −λ j ∆tn c n,j + 1 where j is the precursor group, and γ x,y 0 denotes the zeroth moment while γ x,y 1 denotes the first moment Here x and y are substituted in place of the different types of responses. If properties of the node are fixed within the time step, the zeroth moments in Eq. (10) are equivalent to the total (static) responses, so they can be scored in the same way as explained previously [6,10], by counting the contributions of individual neutrons. The first moments can be scored in a similar way, by multiplying each contribution with the response time via the "internal" neutron clock. Moreover, scoring of the current-to-source, source-to-source, and source-to-current moments needs to be separated according to the type of neutrons contributing to each response (prompt or delayed), similarly to the approach applied in scoring different fission matrices in the TFM method.

TEST CASE
The test model, shown in Figure 1 is based on the Monte Carlo performance benchmark of a Pressurised Water Reactor (PWR) [11]. We have adopted the 17 × 17 fuel pin assembly geometry from the benchmark model and composed a mini-core, consisting of nine such assemblies. The center assembly contains a bank of absorber rods, which can be moved up and down to change the state of the core. The nine assemblies are radially surrounded by a water reflector layer which is one assembly thick. Vacuum boundary conditions are applied on all external surfaces of the model.
We have superimposed the simulation geometry with a mesh consisting of 36 nodes in the z (axial) direction for scoring the response moments and the fission matrices. The responses and the fission matrices were scored for the critical and the delayed super-critical states of the mini-core during separate Monte Carlo criticality calculations. The two states, summarized in Table 1, were obtained by axially moving the absorber bank.
We have simulated an absorber withdrawal transient using TDRM, TFM and TDMC methods. All methods were implemented in an in-house continuous energy Monte Carlo solver. During the first 2 s of the simulation, the set of responses and the fission matrices calculated for the critical state of the mini-core were used; then, at t = 2 s the set of responses and the fission matrices were step-wise changed to the ones calculated for the super-critical state. At t = 6 s the responses and the fission matrices were again step-wise changed to the critical set. In the TDMC simulation, the system was in the critical state during the first 2 s of the transient, then at t = 2 s the tip of the   33.9 ± 0.1 µs 33.8 ± 0.1 µs absorber bank was step-wise moved to the super-critical position. At t = 6 s the absorber was returned to the critical position. The TDMC results were obtained by tracking 4 · 10 5 neutrons, using a 0.2 ms time step for population control. The transient was simulated until t = 10 s.

RESULTS AND DISCUSSION
The evolution of the relative core power during the absorber bank withdrawal transient is shown in Figure 2. The results of the TDRM solution are superimposed over the TFM and the TDMC solutions. All solutions display a prompt jump immediately after the perturbation was applied, followed by an exponential growth. At t = 6 s a negative prompt jump is seen, corresponding to the system being returned to the critical state, followed by an asymptotic approach to the new  Figure 3 shows the TDRM and TFM solutions of the relative power distribution along the axial dimension of the mini-core at several time instances during the simulation. Initially (curve corresponding to t = 1.00 s) the shape of the curve displays a bottom-peaked profile with the maximum at approximately 140 cm from the bottom of the core. Following the perturbation, at t = 6.00 s, the distribution shows a profile with the maximum at approximately 200 cm. Then, following the return to the critical state, the profile returns to the initial shape, and continues to reduce in amplitude.

CONCLUSIONS
In this paper we compared the time-dependent response matrix solution with the transient fission matrix and the time-dependent Monte Carlo solutions for a control rod movement problem in a mini-core reactor geometry. The three solutions were consistent with each other; the transient fission matrix solution displayed a slight under-shoot of the power after the system was returned to the critical state, following a super-critical excursion. The reason for the under-shoot is to be identified. Moreover, further investigations of the time-dependent response matrix method are necessary, including analyses of error and rigorous benchmarking. We plan to implement the method into an established Monte Carlo reactor physics code for this purpose.