SP3-coupled global variance reduction method based on RMC code

Qing-Quan Pan ?Teng-Fei Zhang ?Xiao-Jing Liu?Hui He?Kan Wang

Abstract A global variance reduction(GVR)method based on the SPN method is proposed.First,the global multi-group cross-sections are obtained by Monte Carlo(MC)global homogenization.Then,the SP3 equation is solved to obtain the global flux distribution.Finally,the global weight windowsare approximated by theglobal flux distribution,and the GVR simulation is performed.This GVR method is implemented as an automatic process in the RMC code.The SP3-coupled GVR method was tested on a modified version of C5G7 benchmark with a thickened water shield.The results show that the SP3-coupled GVR method can improve the efficiency of MC criticality calculation.

Keywords RMC code.Global homogenization.Variance reduction.SPN theory

1 Introduction

The reactor criticality calculation is an essential step in thephysical design of anuclear reactor.Thereactor burnup and neutronic/thermal(N–TH)coupling calculations[1]are based on the reactor criticality calculation.If the accuracy of criticality calculation cannot beguaranteed,the accuracy of the reactor burnup and N–TH coupling calculations will also not be guaranteed.The Monte Carlo(MC)method is widely used in reactor criticality calculations because of its high fidelity.However,a non-uniform distribution of relative errors has been found in MC criticality calculations.This non-uniform distribution of relative errors has a negative influence on the efficiency and accuracy of somelocal tallies.In somecases,it even causes numerical instability in MC multi-physics simulation[2]and MC kinetic calculation[3].To improve the efficiency and accuracy of MC criticality calculations and ensure the numerical stability of MC multi-physics simulations,the global variance reduction(GVR)methods should be studied.

There are many variance reduction(VR)methods,such as Zheng’smethod[4,5],theadaptivemethod[6],the FWCADISmethod[7],theautomated weight windowsmethod[8],and the importance function method[9].These methods are well-known for their value in VR applications.However,these methods were initially developed for MC fixed-source calculations and not for MC criticality calculations.Some VR methods have also been developed for MC criticality calculations.They include the UFS[10],UTD[11],and UVM[12]methods.These methods are easily implemented and have good universality.However,the UFS method assumes that the fission neutrons are uniformly distributed in the regions with fuel and have little influence on the non-fuel regions.Moreover,only the sourcebias technique and not theweight window technique is adopted in the UFS method.This limits its value for GVR.The UTD method can reduce the variance of various statistics,but because itsbasic principle is still the same as that of the UFSmethod,it faces the same limitations.The UVM method is newly proposed and can improve the efficiency of MC criticality calculation and flatten the global variance.However,the development of this method is still in the initial stage and needs to be further studied.

A SP3-coupled GVR method is proposed here to improve the accuracy and efficiency of MC criticality calculation.The SP3-coupled method can beused for GVR in MC criticality calculation and allows the analysis of the reactor from its interior to its exterior through a single MC criticality calculation.

The SP3 method is a numerical method for solving the neutron transport equation with a simplified treatment of spatial anisotropy and has been widely studied[13,14].The traditional SP3 method starts from a one-dimensional model.The one-dimensional governing equation is first found and the one-dimensional operator is subsequently replaced with a three-dimensional operator to obtain the three-dimensional governing equation.Compared with the traditional transport theory,the spatial anisotropy is simplified in the SP3 method to achieve higher calculation efficiency.Compared with the traditional diffusion theory,the SP3 method provides a better treatment of spatial anisotropy to achieve higher calculation accuracy.Therefore,the SP3 method balances accuracy and efficiency and has great potential for coupling with the MCmethod to achieve accelerated calculations and GVR.

In this paper,a hybrid MC method based on the SP3 theory is studied and a GVR method driven by the SP3 theory is developed for the RMC code[15].Unlike other hybrid MC methods,such as the fission matrix[16]and response matrix[17]methods,the neutron transport equation is solved instead of a matrix relationship.The remainder of this paper is organized as follows:Sect.2 introducesthetheory and providestherelated formulasand details.Section 3 presents the numerical results with some basic analysis.A discussion is given in Sect.4 before the conclusion.

2 Methods

The SP3-coupled GVR method is implemented as an automatic processin RMCcode.Thisimplies that only the input card of the RMC code is required,and the GVR calculation will subsequently be automatically completed by RMC code.The flowchart for the calculation is shown in Fig.1.The calculation consists of three main steps:

(1) The global cross-sections are determined by MC global homogenization.

(2) The SP3 equation is solved to obtain the global flux distribution.

(3) The global weight windows are determined by the global flux.

Each of the three steps is described in the following subsections.

2.1 MC global homogenization

The SP3-coupled GVR method is a hybrid MC method in which the MC method is coupled with the SP3 method.A deterministic reactor code[18]for solving the SP3 equation,NLSP3,is developed.In contrast to the discrete ordinate method(SN),global multi-group cross-sections(MGCS)are required to solve the SP3 equation.Therefore,the MGCS must first be determined.MC global homogenization was adopted to determine the MGCS for the following reasons:

(1) Because a continuous-energy nuclear library is used,no resonance effect treatment[19]is required.

(2) Because the homogenization process is not based on an assembly model,no leakage correction[20]is required.

(3) Because the geometric model is not approximated,a higher accuracy can be achieved compared with the traditional assembly homogenization.

(4) Because the global homogenization is based on the entire reactor,the cross-sections of all the material regions can be determined through a single calculation.

In particular,Point(4)implies that MC global homogenization is suitable for an automated process.Because the MGCShave a significant influence on the global flux distribution,a high accuracy of the MGCS must be ensured.Meanwhile,because the cross-sections in all the material regions are determined simultaneously,a non-uniform distribution of relative errors will lead to different crosssection precisions in different material regions.To improve the accuracy of cross-sections and ensure that the crosssectionsin all material regionshavethesameprecision,the following measures are taken[21]:

(1) A fixed-source calculation with a flat neutron source is performed for MC global homogenization.The spatial distribution of neutrons is positively correlated with the neutron flux distribution if the VR techniqueisadopted in the MCsimulation.In alarge reactor,the neutron fluxes in different phase spaces may differ,especially when the reflective and shielding layers are considered.These differences result in a variation in the spatial distribution of neutrons.Therefore,although performing global homogenization is a reasonable approach in MC criticality calculation,it is inevitable that the crosssections in some regions will have higher precision compared to other regions.

The spatial distribution of neutronscan be artificially set in a fixed-source calculation.Suppose a flat neutron source is provided in the fixed-source calculation.In this case,all the material regions have almost the same density of neutron histories,and the cross-sections obtained will all have almost the same precision.Therefore,the MCglobal homogenization in the fixed-source calculation is performed with a flat neutron source.

(2) The neutron fission reaction is restricted during neutron transport.Various reactions such as scattering,fission,and absorption occur during neutron transport.The reactor becomes critical when the loss of a neutron because of the absorption of neutron or its transportation outside the system is compensated by the corresponding production of a neutron in a fission reaction.Typically,a fission reaction produces two or three neutrons.

In a supercritical system,the loss of one neutron is accompanied by the production of more than one neutron.Therefore,if fission neutrons are produced during the fixed-source calculation,the number of simulated neutrons will increase,and the calculation will never be terminated.As stated earlier,global homogenization is performed during the MC fixedsource calculation.To ensure the feasibility of the fixed-source calculation,neutron fission reaction is restricted during MC global homogenization:The neutrons only perform random walks in the system,and no fission reaction occurs.However,various reaction probabilities are recorded during neutron transport.

(3) The critical spectrum isapproximated using the Watt fission spectrum.The MC homogenization involves a process of compressing and grouping the neutron energy spectrum.The neutron spectrum is essential in MCglobal homogenization.However,becausethe critical spectrum isunknown before the fixed-source calculation,it can only be approximated.

As stated before,neutron fission is restricted during the neutron transport process.This restriction leads to errors.Therefore,the neutron spectrum can be artificially adjusted according to the fission spectrum when the flat neutron source is set during the fixedsource calculation.It is known that the fission spectrum of U-235 is described by the Watt spectrum[22].Therefore,the external fission source is distributed according to the Watt spectrum.

where C,a,and b are coefficients that are ordinarily equal to 0.453 MeV,0.965 MeV,and 2.29 MeV-1,respectively.

It should be noted that although the Watt spectrum may not be a good one for the external source and is not universal for different fission nuclides,it is the best approximation that we can come up with now.

(4) Corrections are adopted to make the neutron spectrum asymptotic to the criticality spectrum.The main consideration in the Watt fission spectrum are the fission reactions that occur during the neutron transport process.The Watt fission spectrum is a good approximation in the reactor core region where many fission reactions occur.However,in regions without fissile materials,the neutron spectrum is mainly influenced by scattering reactions,and the error will be larger when the Watt spectrum is used to approximate the critical spectrum.Therefore,there should be as many scattering reactions as possible before the neutrons are lost.After many collisions,the neutron energy spectrum in the reflective and shielding layers become close to the real spectrum.To make the neutron spectrum asymptotically approach the critical spectrum,the cutoffweight limit is reduced when more scattering reactions occur.In the SP3-coupled GVR method in this study,the weight cutoffthreshold was set to 1.0E-8,which isthedefault value in the RMCcode for the VRsimulation.In thenonanalog MCmethod,the scattering reactions are effectively increased by reducing the cutoffweight windows.After multiple collisions in the reflective and shielding layers,the energy spectrum of the neutrons becomes close to the critical spectrum.However,at this stage,the weight of the neutrons is small.The neutrons therefore have little influence on the tallied results.To increase the influence of neutrons after multiple collisions on the cross-sections,the reaction rate is recorded based on the number of reactions rather than the weight of neutrons.

We calculated the MGCSusing a fixed-source calculation,which is an approximation.Thus,the MGCSobtained inherently has errors,which implies that the efficiency is improved at the cost of calculation accuracy.The SP3 solution only provides the parameters for VR,and the final result is given by MC criticality calculation.Thus,even if the accuracy of MGCS is compromised,the error in the MGCS does not greatly affect the accuracy of the final result.This is one of the reasons why we adopt an approximate method to solve for the MGCS.

The SP3 equation is solved to obtain the global flux distribution.The total grouped cross-sections Σt，g,grouped absorption cross-sections Σa，g,grouped fission production cross-sections νΣf，g,inter-group scattering matrix Σs，g′→g,and fission spectrumχgneed to be input to the NLSP3 code.Therefore,only these five grouped cross-sections are calculated in the MC global homogenization.It should be stated that it is not easy to obtain the anisotropic crosssections during the MC global homogenization.To reduce the computational cost,only isotropic cross-sections are considered.

The grouped flux distribution and grouped reaction rates are tallied to calculate the grouped cross-sections.The traditional volume-weighted method is employed in which the MGCSis defined as

where r is a spatial variable,φg（r）is the global flux distribution,and the subscript i in Vilabels the different spatial regions.

The track length is used to tally the flux distribution and reaction rates.The flux is calculated as

The homogenized cross-section is therefore

where Egand Eg-1are the lower and upper bounds of the gth energy interval,respectively;W.Tis the product of the weight and track length of the ith neutron in region V;and N is the total number of simulated neutrons.

The calculation of fission energy spectrum is relatively simple.The fission energy spectrum is calculated with reference to the number of fission neutrons produced in each energy group.The number of fission neutrons produced is then normalized to obtain the energy-dependent probability of neutron production in the energy interval E.The fission energy spectrum is equal to the energy-dependent probability of neutron production,that is,

where p（E）is the energy-dependent probability of neutron production.

The process of calculating the scattering matrix is the same as that for calculating the fission energy spectrum.The neutron energies before and after scattering are tallied to calculate the scattering probability:

Using the scattering probability,the inter-group transfer cross-section is calculated as

The MGCS is calculated using Eqs.(3)–(8)during the MC global homogenization.The MGCS of all material regions required by the NLSP3 code can thus be determined simultaneously in a single fixed-source calculation for subsequent use in the SP3 global calculation.

2.2 Solution of the SP3 equation

The SP3 equation is solved to obtain the global flux distribution.The SP3 equation is adopted because of the following reasons:

(1) Compared with the diffusion and CMFD equations,the SP3 equation has higher accuracy[23].

(2) Solving the SP3 equation takes approximately the same amount of time as solving the diffusion equation[24].

(3) The SP3 equation is easier to be applied to larger systems because it requires less memory and computing time compared to the SN equation.

(4) To date,there has been no research on coupling the SP3 theory and the MC method for GVR.For academic interest,we would like to make an attempt to investigate this topic.

The SPN theory[13,14]is a simplified spherical harmonics method for solving the neutron transport equation.It is widely studied because it has great potential for nextgeneration reactor physics calculations.There are two variants of the SPN theory,namely,the traditional and the rigorous SPN theories.In the traditional SPN theory,the one-dimensional governing equation is first found from a one-dimensional model,and the one-dimensional operator is subsequently replaced with a three-dimensional operator to obtain the three-dimensional governing equation.The traditional SPN theory lacksarigorousmathematical basis.The rigorous SPN theory is derived based on the variational method and theconcept of grouped angular flux.The angular flux is clearly defined in the rigorous SPN theory so that the physical fundamentals are completely considered in the treatment of boundary conditions.These two SPN theories share the same SPN equation but use different boundary conditions.

Although a complete mathematical derivation is provided in therigorous SPN,thesolution of therigorous SPN equation suffers from numerical instability because of its incompatibility with the traditional transverse integral technique[25].The traditional SPN equation is therefore solved in the SP3-coupled GVR method.The ‘N’in‘SPN’represents the order of the polynomial expansion.When N is equal to 3,the SP3 equations are obtained as:

where the superscript k isthenodal index,the subscript g is the energy interval index,is the zeroth-order scalar flux at the k th node in the gth energy interval,is the second-order scalar flux at the k th node in the gth energy interval,and

The source is

The boundary conditions[26]for the traditional SP3 equation are

where the superscript M represents the quantities at the node boundary,andandare the zeroth-order and second-order partial currents,respectively.

The numerical instability problem was encountered when the SP3 equation was solved using the traditional nonlinear iterative method[27].The zeroth-order flux is usually hundreds of timeslarger than thesecond-order flux.When the traditional coupling corrective relationship was used to calculate the corrective factors,there was also a 100-fold difference between the zeroth-order and secondorder corrective factors,which led to an ill-conditioned iteration matrix and convergence failure.The following new coupling corrective relationship is proposed to solve the numerical instability problem:

A new nonlinear iterative method[25]was proposed based on the above mentioned coupling corrective relationship,and the NLSP3 reactor core code was developed.Figure 2 illustrates the nonlinear iterative strategy.Using the NLSP3 code and global cross-sections,the global flux distribution can be obtained easily.

2.3 GVR for MC criticality calculation

The MC criticality calculation is composed of several inactive and active cycles.Each cycle is a fixed-source calculation with a cycle-estimated result.All the active cycles are averaged as follows to obtain the final result of MC criticality calculation:

where the subscript n denotes the cycle index,the final result of MC criticality calculation,N the total number of active cycles,andthe cycle-estimated result of the nth active cycle.is given by

The standard deviation of the final result is calculated from the deviation between the individual active cycles:

Therefore,the standard deviation increases with the deviation between the individual active cycles.The deviation between the individual active cycles should therefore be reduced to reduce the standard deviation of MC criticality calculation.As stated before,a fixed-source calculation is performed for each cycle,and the deviation between the individual cycles is reduced by various VR techniques.

Fig.2 Nonlinear iterative strategy of the NLSP3 code

Supposethat all theactivecyclesareperformed with the same weight windows.The individual cycle-estimated results will then be close to another.This leads to a reduction in the standard deviation.Therefore,the key to improve MC criticality calculation is to obtain high-precision weight window parameters efficiently.

In the SP3-coupled GVR method,only the forward SP3 calculation is performed.Because there is no known adjoint theory of the SP3 method,we cannot obtain the adjoint flux for the time being.We have therefore focused on GVR for the criticality calculations in the meantime.We believe that the forward flux distribution is sufficient for determining the weight window parameters.Therefore,we use the forward flux distribution to approximate the weight window parameters.

The global flux distribution determined by the SP3 equation is denoted as,where the superscript k denotes the spatial region and the subscript g denotes the energy interval.The SP3-coupled method can only bias the spatial position of neutrons.Therefore,the independent energy variable should be eliminated by accumulating the energydependent flux in the spatial region:

where φkis the space-dependent flux distribution,which is normalized before it is used to approximate the global weight windows,and

where max（φk）is the largest space-dependent flux.Using the normalized global flux distribution,the global weight windows are determined as

where Ws,kis the global space-dependent survival weight window in the k th space region.Using Eq.(20),the weight windows of all spatial regions can be determined,following which the ceiling and floor weight windows can be determined using the following formulas:

and

where Wc,kand Wf,kare the ceiling and floor weight windows in the kth space region,respectively.The standard weight window technique with the weight windows of all spatial regions is then used for GVR.

3 Numerical results

The SP3-coupled GVR method was tested on the C5G7 benchmark[28].In the original C5G7 model,GVR can be easily achieved because the water shield layer is not sufficiently thick to pose a deep penetration problem.To better demonstrate the advantages of the SP3-coupled method for GVR,we artificially increased the thickness of the shield layer.The materials and geometry of the reactor core were derived from the C5G7 benchmark,and the core region was surrounded by a water shield layer,which was 64.362 cm thick and three timesthicker than the side of the core assembly side as shown in Fig.3.

Fig.3 (Color online)Modified C5G7 benchmark

Thesuperimposed mesh geometry wasused for the SP3-coupled GVR simulation.The C5G7 reactor core had a 4×4 assembly layout and the assemblies was surrounded by a water shield that was three times thicker as the assembly.The model was divided into a total of 10×10=100 material meshes in the modified C5G7 benchmark.The MGCS of the 100 meshes were determined through the MC global homogenization over the four energy intervals of[0,6.25 × 10–7],[6.25 × 10–7,5.53 × 10–3],[5.53 × 10–3,8.21 × 10–1],[8.21 × 10–1,[20]MeV[29].

A fixed-source calculation was performed for the MC global homogenization.The choiceof the number of source neutrons per material mesh is a challenging issue.A larger number of source neutrons for the fixed-source calculation will increase the precision of obtained cross-sections but will also require more computation.To achieve a balance between efficiency and accuracy,the number of source neutrons per material mesh was chosen as 10,000 based on our experience.

Therefore,a total of 1,000,000 source neutrons were simulated for the MC global homogenization.The SP3 solution only provided weight windows,and the final result was given by MCcriticality calculation.Thus,even if there was an error in the MGCS,the accuracy of final result was not affected greatly.This is why the number of source neutronsper material mesh wasempirically determined for solving the MGCS.

The RMC code called the SP3 solver for global calculation with the cross-sections of the 10×10 material meshes.Each material mesh was divided into a total of 4×4=16 nodes with each node being 5 cm×5 cm in size.Therefore,atotal of 40×40=1600 nodeswere used to solve the SP3 equation.After solving the SP3 equation,the flux distribution in the 1600 nodes was obtained.The three convergence criteria for the inner,outer,and keffiterations in the NLSP3 code were set to 10–3,10–3,and 10–4,respectively.The calculations were performed with 24 threads,and the calculation time taken for MC global homogenization and SP3 solver was 1.14 mins.

Given the global flux distribution,the global weight windows can be determined using Eqs.(20)–(22).The MC criticality calculation was then performed using the global weight windows.There were 10,000 neutronsper cycle,50 inactive cycles,and 450 active cycles in the SP3-coupled GVR criticality calculation.15.67 mins were required to complete the MC global homogenization,SP3 solver solution,and GVR criticality calculation.

For comparison,the standard MC simulation method was also performed to calculate the flux distribution of the 1600 nodes.Because GVR was not performed in the standard simulation method,we call this the direct simulation.Therewere 1000,000 neutronsper cycle,50 inactive cycles,and 450 active cycles in the direct simulation.The calculation time was 151.80 mins.The distributions of relative errors in the SP3-coupled GVR and the direct simulations are compared in Fig.4.The distributions of total flux in the SP3-coupled GVR and the direct simulations are compared in Fig.5.

As shown in Fig.5,the distributions of total flux obtained by the direct and the SP3-coupled GVR simulations are in good agreement.The deviation was calculated to be less than three times the standard deviation.Meanwhile,it can be seen from Fig.4 that when the modified C5G7 benchmark wassimulated directly,therelativeerrors increased with the distance away from the reactor core.Compared with the direct simulation,the global variance distribution in the SP3-coupled GVR simulation was flattened in the entire system.High accuracy was achieved in both thereactor coreand theshielding layer.Therefore,the SP3-coupled GVR method can improve the efficiency of MC criticality calculations.To better quantify its GVR performance,three quantities were compared:

Fig.4 (Color online)Variance distribution of the SP3-coupled GVR and direct simulations

Fig.5 (Color online)Distribution of the total flux in the SP3-coupled GVR and direct simulations

1.The average figure of merit(AV.FOM),which is defined as

2. The standard deviation of the relative deviation,which is defined as

3.The ratio of the maximum to the minimum variances(σmax/σmin),i.e.,

The quantitative results for the direct and SP3-coupled GVR simulations are compared in Table 1.

Table 1 Results of GVR and direct simulations

The results show that,compared with the direct simulation,the average figure of merit(AV.FOM)was improved 33 folds in the SP3-coupled GVR method.Meanwhile,the ratio of the maximum to minimum variances was reduced by a factor of five.The value of the SP3-coupled GVR method for GVR and its ability to improve the efficiency of MC criticality calculation are hence clearly demonstrated.

It has been argued that it is unnecessary to couple the SP3 method with the MC method for GVR,and that the diffusion method is sufficient.To demonstrate the advantage of coupling the SP3 method with the MC method for GVR over coupling the diffusion method,we also developed a diffusion-coupled GVR method based on the RMC code.The flowcharts of the diffusion-coupled GVR and the SP3-coupled GVR methods are almost identical except for the solver used to solve the global flux distribution.In the SP3-coupled GVR method,the SP3 equation is solved,while in the diffusion-coupled method,the diffusion equation is solved.To solve the diffusion equation,the diffusion coefficient is required.Therefore,the total cross-section Σt，gused in the SP3 equation was converted into the diffusion coefficient Dg.

Calculations were then performed on the same model using the diffusion-coupled GVR method with the same calculation parameters.The obtained AV.FOM from the results of the diffusion-coupled method was 136.81.As shown in Table 1,the AV.FOM of the result of the SP3-coupled method was 383.94,which is2.81 timesthat of the diffusion-coupled method.Therefore,the SP3-coupled GVR method showed better GVR performance than the diffusion-coupled GVR method.We believe that the SP3 method has great potential to be coupled with the MC method for accelerated calculations and GVR.

4 Discussion and conclusion

An automatic process based on the MC global homogenization,the SP3 global solver,and GVR criticality calculation was developed in the SP3-coupled GVR method.The numerical results show that this method can improve the efficiency of the MC criticality calculation.

The SP3-coupled method proposed in this paper was developed for MC criticality calculations.It,therefore,superficially resembles the traditional two-step reactor physical calculation framework.However,the SP3-coupled method differs from the traditional framework in three aspects:(1)global homogenization isperformed,instead of assembly homogenization;(2)the homogenization process is completed during the fixed-source calculation,and not during the criticality calculation;and(3)global MC homogenization and the global reactor calculation are achieved in a one-step process instead of a two-step process.Moreover,only the weight window parameters for MC criticality calculation are output in the SP3-coupled method.The final results are obtained from the MC calculation rather than the global calculation of the NLSP3 code.

The SP3-coupled method has similarities and differences with other hybrid MC methods,such as the fission matrix and the response matrix methods.They share similar calculation processes in which some local information of the system is recorded in advance and used to build a global formula,and the global information is obtained by solving the global formula.In the fission matrix method,the fission coupling relationship between different regions is recorded first.In the response matrix method,the relationship between the currents and sources in different regions is recorded first.In the SP3-coupled method,the homogenized cross-sections in different regions are recorded first.The matrix relationship is solved in the fission matrix and response matrix methods to obtain global information.In contrast,in the SP3-coupled method,the neutron transport equation issolved.We,therefore,believe that the SP3-coupled method has advantages in terms of numerical stability,development potential,and application scope.

In this study,the traditional SP3 equation is solved in the SP3-coupled method and the forward flux is used to approximate the weight windows.To further improve its performance for GVR,the rigorous SPN theory proposed by Prof.Yung-An Chao and the adjoint method for the SPN theory should be studied to obtain better weight window parameters.

Author contributionsAll authors contributed to the study conception and design.Material preparation,data collection,and analysis were performed by Qing-Quan Pan,Xiao-Jing Liu,and Kan Wang.The first draft of the manuscript was written by Qing-Quan Pan and all authors commented on previous versions of the manuscript.All authors read and approved the final manuscript.