A Hybrid Monte-Carlo-Deterministic Method for AP1000 Ex-core Detector Response Simulation

The ex-core detector-response calculation is a typical deep-penetration problem, which is challenging for the Monte Carlo method. The response of the ex-core detector is an important parameter for the safe operation of the nuclear power plants. Meanwhile, evaluation of the ex-core detector response during each step of fuel-loading is used to guide the fuel-loading sequence. The response can also be used to reconstruct corepower distribution for online monitoring of long-term power. The detector used for the ex-core response is the source-range detector which is sensitive to thermal neutrons. For a Monte Carlo shielding calculation of the above detector response, the thermal flux under 0.625eV is needed, which is too small to be tallied by traditional Monte Carlo simulations. In practice, the tally results are close to zero in the detector region under direct Monte Carlo calculation. Even if the number of particles is increased to a significant amount, the statistical variance is still very large. The high variance along with a significant calculation time leads to a small Figure Of Merit (FOM). In order to solve this problem and to improve the tally efficiency of the ex-core detector response, a hybrid Monte-Carlo-deterministic method is employed in this study, and an in-house hybrid Monte-Carlo-deterministic particle transport code, NECP-MCX, is developed in this paper. The method takes the space-energy-dependent adjoint fluxes to generate importance parameters for the mesh-based weight window in the Monte Carlo calculation. Simultaneously, the mesh-based source biasing is performed with the consistent importance parameters to make the starting weight of neutrons matching with the survival weight of the weight windows. As the mesh used in the hybrid MonteCarlo-deterministic method is superimposed, the mesh of the weight window will not be affected by the complex geometry model. The adjoint flux is obtained by the efficient SN method with the multi-group cross-section data. The whole toolset is convenient to use with single set of the modelling data for both Monte Carlo and deterministic simulations. Compared with the direct Monte Carlo simulation, the hybrid Monte-Carlo-deterministic method has a higher efficiency for a typical deeppenetration problem such as the AP1000 ex-core detector-response simulation.

AP1000 is a 1000MWe pressurized water reactor (PWR) with passive safety features. During its commissioning/operation time, the response of ex-core detector is important for getting the information of commissioning/operation state. The ex-core detectors of AP1000 are located outside the thermal insulation layer which are far away from the reactor core. There are thermal insulation layer, pressure vessel, neutron pad, core barrel and core baffle between the detectors and the core, which have great effect of shielding. As a result, the flux will decrease about 10 orders of magnitude, making the detector signal very weak if the source is not strong enough. On the other hand, the ex-core detectors of this reactor are source-range detectors which are sensitive to the thermal neutrons only [1,2], which is more difficult for the detectors to get the signal. However, an accurate calculation of the detector response is essential for the core design and the state prediction before the commissioning/operation of the reactor [3].
The arrangement of the inner core of AP1000 is complex [3]. The 18-month-cycle low-leakage design uses five fuel regions with U 235 enrichment ranges from 0.7 w/o to 4.8 w/o. Enrichment zoning and burnable absorbers are used for power-shape and reactivity control. In particular, the use of the Integral Fuel Burnable Absorber (IFBA), a thin ZrB 2 coating on the pellet surface, and the discrete Wet Annular Burnable Absorber (WABA), inserted at selected guide thimble location, makes the detail AP-1000 whole-core modelling more complex. Great challenges have been encountered when using traditional deterministic shieldingcalculation codes such as discrete-ordinates codes [4,5] to model the AP1000 whole-core with such complex structure to get an accurate ex-core detector response.
On the other hand, the Monte Carlo method has great advantages over the deterministic method in shielding calculation especially for the modelling of complex geometry structure such as AP1000. During the commissioning of AP1000, two small start-up sources made of californium were set inside the core. The start-up sources, with a height of 3.8 cm, are put in the control-rod guide tube of one fuel assembly. They are equivalent to point sources relative to the scale of the whole core. The ex-core detectors are located far away from these two sources. With all assemblies in the core, a particle from the source should go through many facilities to reach the ex-core detectors. It is a typical deep-penetration case with low tally or even zero tally in the target region. It is hard to directly simulate particles and get results with relatively low statistical error with a Monte Carlo code. In the past, many variance-reduction methods such as the region importance, the weight window, cutoff of energy or weight, and the source biasing, etc. have been developed to handle this kind of problem. The usage of these methods is empirical which needs extensive experiences and numerical iterations in order to get acceptable results. They are not convenient to use in lack of experience (e.g. the simulation of AP-1000 ex-core detector response during the commissioning phase). Sometimes the efficiency may get worse if biased in a wrong direction. A method was proposed in reference [6] which uses the adjoint-flux distribution as the importance map of the phase space for the Monte Carlo simulation.
In this paper, a hybrid method with mesh-based weight window is employed. A new hybrid Monte-Carlodeterministic particle transport code is developed at the NECP lab in Xi'an Jiaotong University. The computer code, called NECP-MCX, embeds the in-house S N code NECP-hydra [5] to perform the adjointflux calculation as part of the hybrid method. The importance theory is based on the consistent adjointdriven importance sampling (CADIS) theory [6,7], in which the adjoint flux is used to set the weight window and to bias the source distribution. NECP-MCX is designed to input only one set of data that contains the traditional CSG modeling information for the Monte Carlo code. Special algorithm has been developed to generate the S N model, and to match mesh bins between the source distribution and the S N calculation. NECP-MCX can set the weight window and the source biasing automatically and it performs well with very high Figure Of Merit (FOM) when compared with the direct Monte Carlo simulation or the Monte Carlo calculation with the region-importance option. With NECP-MCX, the AP1000 ex-core detector response can be simulated accurately and efficiently without users' experiences.

THEORY OF THE HYBRID MONTE-CARLO-DETERMINISTIC METHOD
In the Monte Carlo calculation of the neutron-transport equation, a particle with the information of location, direction, energy and weight is simulated as it is in the real physical world. It goes through the material and collide with the target nuclide. While for deep-penetration cases, with strong absorbent material among the case, particles are difficult to reach the region far away from the source. As a result, the sample space which is the particle number in the target region is very small making the variance of the tally result too large. More source particles are needed to increase the number of particles in the target region. The increasing of source particles will increase the simulation time proportionally. Under the theory of probability, when the variance of tally result increases an order of magnitude, the number of particles should increase two orders of magnitude, so does the simulation time. Thus, simply increase the particle number is not an efficient way to reduce the variance of the tally result. And for the evaluation of the efficiency, the Figure Of Merit (FOM) [8], as shown in equation (1), is used: R is the statistical relative error of the tally, which is the standard deviation divided by the average value. T is the total calculation time. Actually, FOM does not change with the increasing of the particle number.
The main idea of reducing the variance of the target region is to make the population of the particle large. Under this criterion, the bias of transport and the bias of source are applied. In the traditional way to set the weight window, a set of weight-window parameters are applied in a certain phase space. These parameters include the upper bound WU , the lower bound WL , and the survival weight WS . When a particle enters into this phase space, the weight of particle wgt is checked. If the weight is higher than the upper bound of the weight window, split the particle into [ ] / wgt WU secondary particles. If the weight is inner the weight window, continue to simulate without special process. If the weight in lower than the lower bound of the weight window, play a Russia Roulette to the particle with the probability of / wgt WS , the survival particle will have a weight of WS . With this logic of weight-window checking, the number of particles will increase in the region with lower weight window, which means the importance is higher in this region.
How to determine the weight-window parameter is a problem. In the consistent adjoint driven importance sampling theory, the adjoint flux is used to set the weight window and bias the source distribution.
According to the properties of the adjoint equation, the forward flux ( , ) r E φ and the source ( , ) q r E has the relation with the adjoint flux * ( , ) r E φ and the adjoint source * ( , ) q r E . The relation is shown in the following equation: The angle brackets "< >" means integrating over the whole phase space of r , E and Ω .
The following equation is used to get the detector response in a fixed-source calculation: Σ is the response function of the target tally which is zero outside the target region. Normally the flux of the target region is calculated and the integral of the flux and the response function in the target region is the detector response RES . Now, if we set the response function d Σ as adjoint source * q , we obtain the following alternative expression for the detector response: The detector response can be obtained by the integral of the forward source q and the adjoint flux * φ .
From the equation (4), it can be concluded that the adjoint flux * ( , ) r E φ represents the contribution of the source to the detector response. Under this consideration, take the distribution of the adjoint flux as the weight-window parameter is the most proper way to define the weight window. When the adjoint flux is higher in one region, the weight window will be lower compared with other regions. Particle that enters this region will split into more population which makes the sample space larger in this region, thus the variance can be reduced.
In order to apply the above weight-window setting and the source biasing, the adjoint flux is needed. However, the solution of the adjoint transport equation is as difficult as the forward transport equation. It seems to be an endless loop to solve this problem. Since the adjoint flux is used to guide the importance distribution, the requirement of the accuracy is not too high. Introducing some approximation in the adjoint transport calculation is acceptable. Hence the S N method is chosen to solve the adjoint transport equation with some modelling approximation such as use of the multi-group nuclear data and coarse-mesh geometry to derive a coarse importance map.
Under the theory above, the adjoint flux is obtained from the S N solver. So, the distribution of the adjoint flux is based on the mesh in space and based on the group in energy. In the real case, the structure of the source distribution and the importance map cannot be exactly same with the structure of the adjoint-flux distribution, neither in space nor in energy. The mapping among the source distribution, the importance distribution and the adjoint-flux distribution is needed.
First is the distribution of the energy, each mesh bin in space has a set of energy distribution with a same group structure. There are three sets of energy-group structures in the hybrid code. First is the energy-group structure of the source, it depends on the source definition. Second is the energy-group structure of the S N solver based on the multi-group nuclear-data library. In this paper for the calculation of the AP1000 ex-core detector, the Bugle96 library with 47 energy groups was used. Third is the energy-group structure of the weight window. The distribution of each energy group is regarded as a flat distribution, so the mapping of the group distribution is under the weight by the energy range.
In the above equation, 1 j e − ∆ is the energy range of ( 1) j th − energy group of the input energy structure.
is part of the energy range of the -i th energy group of the output energy structure that overlaps in the j-th energy group of the input energy structure. With this kind of algorithm, the distribution of the output energy structure can be got from the distribution of the input energy range.
Similar to the energy-group structures, there are also three sets of spatial mesh structures in the hybrid method for the definition of the space distribution. All meshes are in the rectangular form. The same algorithm used for the energy-group mapping is applied in the spatial-mesh mapping. Because the spatial mesh is in three dimensions, three nested structures are used to map the input spatial-mesh structure to the output spatial-mesh structure.
Based on the method and the algorithm described above, a hybrid Monte-Carlo-deterministic particletransport code NECP-MCX has been developed. NECP-MCX is designed as a compact toolset which minimizes users' effort. The flow of the hybrid Monte-Carlo-deterministic method in NECP-MCX is show in the Fig. 1.

DEVELOPMENT OF THE HYBRID MONTE-CARLO-DETERMINISTIC CODE
To address the challenges of deep-penetration calculations faced by the conventional MC codes, a new hybrid Monte-Carlo-deterministic particle-transport code NECP-MCX has been developed. NECP-MCX combines the advantages of the MC method and the S N method. The code utilizes an in-house S N code NECP-Hydra to quickly perform the forward and adjoint transport calculation with higher efficiency. The adjoint flux obtained is then used to set consistent weight-window and source-biasing parameters which are employed to reduce the variance in Monte Carlo (MC) simulations for deep-penetration problems.
NECP-Hydra is a three-dimensional parallel discrete ordinate transport code developed for high parallel efficiency and scaling efficiency on large number of processors. Many advanced and useful features are included in NECP-Hydra to make the code more practical in the engineering application. NECP-Hydra has reliable accuracy for both the fixed source and eigenvalue problems.
NECP-MCX is able to simulate neutron and photon transport based on either the continuous-energy library or the multi-group library. NECP-MCX has the functionality to simulate both the fixed-source problem and the eigenvalue problem. NECP-MCX applies classic variance-reduction techniques such as region importance, weight window and source biasing. To obtain reliable values of the tally parameters and associated variance, NECP-MCX simulates several batches with many particle histories for each batch.
Except the capability of deep-penetration simulation, the capabilities of neutron-photon coupling calculation and criticality calculation based on the Monte Carlo method have also been developed in NECP-MCX. The capabilities of burnup, activation, source-term and multi-physics calculation are currently under development for wider applications of the code in the future.
The state-of-art acceleration techniques for the MC method such as the hash-based energy-lookup algorithm, the hash-based mapping algorithm, the nearest-neighbors algorithm, etc. have also been implemented in NECP-MCX. The whole toolset is convenient to use with single set of the modelling data for both Monte Carlo and deterministic simulations.

APPLICATION IN AP1000 EX-CORE DETECTOR CALCULATION
AP1000 is a large pressurized water reactor designed by Westinghouse. There are 12 kinds of different assemblies in the core. They are under the low-leakage design with five fuel regions. In this paper, the response of the ex-core detector were calculated with all fuel assemblies loaded into the core. Before the criticality of the reactor, the start-up sources provide primary neutrons to induce fission in the core. After loading all fuel assemblies into the core, two Primary neutron-source assemblies are in the core to provide neutrons. The arrangement of these two assemblies are shown in Fig. 2. The location of the source pin is marked as black in Fig. 2. The source is made of californium, with the height of 3.8cm. It is 106.68cm away from top boundary of the bottom reflector. The source pin is consist of He-4 gas, Cf-252 source, the SS-304, the boron water and Zr cladding. In the calculation reported in this paper, the source was defined as a pin box with the intensity of the source intact. Considering the scale of whole core, this approximation is acceptable. The location of primary neutron-source assemblies is shown in Fig. 3 (a).
Outside the fuel assembly, there are core baffle, core barrel, neutron pad, pressure vessel and thermal insulation layer. The structure is shown in Fig. 3(b). The detectors are outside of the thermal insulation layer. There are totally four detectors in the south, north, west and east of the core respectively, which named as A, B, C and D in the analysis of this paper, as is shown in Fig.3(b). All of them are 3.213 m away from the center of the core. The detector was modelled as a cylinder region with a cladding. The flux under 0.625eV was tallied to get signal of detector. Because the tally region of the detector is far away from the source which makes it difficult to tally, especially for the flux under 0.625eV. It is a typical deep-penetration problem with very high variance even with the variance-reduction technique.
To show the efficiency of the hybrid Monte-Carlo-deterministic method, the NECP-MCX results are compared with those of MCNP with the imp option [8]. Imp in MCNP means region importance set by users manually based on users' experience for different cases. The measured results from nuclear plant is used as the reference. The results are in the form of cps, which means count per second. It is worth noting that since the count at the detector is small, the measured results fluctuate greatly.
In the simulation of this case, totally 2 billion particles were simulated to get a reliable tally result in detector region. The nuclear-data library used by NECP-MCX is the ENDF-B-VII.0 continuous-energy library. The nuclear-data library used by NECP-hydra is the BUGLE96 multi-group library for shielding calculation.
The mesh structure of S N calculation is shown in Fig. 3(c). Finer mesh was used in the detector regions in order to get more accurate adjoint-flux distribution around detectors. The S N order is S 8 in the simulation. With normal input of the Monte Carlo model and with the hybrid switch turned on, the code will automatically generate the model of S N solver and do the adjoint calculation first. Then the adjoint flux of each mesh and each group was used to set the weight window and to bias the source after mapping of groups and meshes as described in Section 2. Finally, the Monte Carlo simulation will be performed under the new weight window and the biased source.

CONCLUSIONS
For the commissioning of AP1000, the two start-up sources are the primary sources in the core. With the full loading of fuel assemblies, the response of ex-core detector is very difficult to calculate with either the deterministic method or the Monte Carlo method. A hybrid method based on the importance theory is applied in this paper, which automatically generates the importance distribution from the S N model and performs the mapping of the distribution for each energy group and each mesh bin. The hybrid Monte-Carlo-deterministic particle-transport code NECP-MCX has been developed and applied in the calculation of AP1000 ex-core detector response. The response results of AP1000 ex-core detectors simulated by NECP-MCX performs very well with small statistic relative errors. Compared with the poor (almost zero) tally results from the direct Monte Carlo simulations, the performance of the hybrid code NECP-MCX is remarkable. Compared with the results from the MCNP calculation with the optimized imp setting based on experiences, the Figure Of Merit (FOM) increases significantly. In addition, since the hybrid code NECP-MCX automatically sets the weight window and biases the source based on the S N simulations, it is much more convenient to use NECP-MCX than MCNP which needs manual imp setting based on experiences. With NECP-MCX, the deep-penetration problem such as AP1000 ex-core detector response can be simulated accurately and efficiently without users' experiences.