Nuclear data

The data were from the ENDF/B-VII ^[17] and CENDL-3.2 ^[18] evaluation libraries. The cross-section parameters were generated by using NJOY code^[11]. The energy ranges are 10^-11-20 MeV for neutrons, 1 eV-100 GeV for photons, and 10 eV-1 GeV for electrons, especial for 1 eV-1 keV of photons and 10 eV-1 keV of electrons are newly increased. It makes that JMCT has a simulation capability for visible light problems. The continuous point-wise cross-section accounts for all the neutron reactions in the evaluation library. Thermal neutron processes in free-gas mode or S(α,β) mode. For photons, coherent and incoherent scattering, fluorescence emission, and electron pair generation are considered to follow photoelectric absorption. Nuclear data tables exist for neutron interactions, neutron-induced photons, photon interactions, neutron dosimetry or activation, and thermal scattering, S(α,β). Photon and electron data are atomic rather than nuclear in nature. The functions of proton, atmosphere and molecule transport are newly added in JMCT3.0.

Simultaneously, JMCT provides the function of multigroup calculation, where a 47-group neutron and 20-group photon P₅ cross-section library is mainly used for radiation shielding, which its energy structure is same as the Bugle library ^[19]. Another 172-group neutrons, 40-group upper scatter neutrons, and 30-group photon P₅ library is mainly used for the reactor core.

Source specification

JMCT provides several standard sources. Users can specify a wide variety of source conditions. Independent probability distributions may be specified for the source variables of the position, energy, direction, time, and zone or surface. Information regarding the geometrical extent of the source can also be provided. In addition, the source variables may depend on other source variables. The user can bias all of the input distributions. For the energy it includes various analytic functions for fission and fusion energy spectra, such as Watt, Maxwellian and Gaussian spectra. Biasing can also be accomplished by using special building functions. Simultaneously, the user source subroutine also provides a standard port.

Estimation of errors

Monte Carlo error e satisfies the center limit theorem as following $\underset{N\to \infty }{\lim }P\left\{\left|\frac{\overline{x}-a}{\overline{x}}\right|<\epsilon \right\}=\Phi (\lambda )$ (1) where $\overline{x}=\frac{1}{N}{\displaystyle {\sum }_{i=1}^N{x}_i}$ is the mean, N is the sample number (or histories), $a=E(x_i),i=1,2,\cdots$ is a mathematical expectation of various random x_i, $\epsilon =\frac{\lambda \sigma }{\sqrt{N}\overline{x}}$ is an error, s σ is a variance and $\Phi (\lambda )=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\lambda }^{\lambda }{e}^{-t^2/2}dt}$ is the degree of confidence.

$\Phi (1)=0.6827$, $\Phi (2)=0.9545$, $\Phi (3)=0.9973$. In general, the error in $\Phi (1)$ is the standard deviation of the fixed-source problem. The errors in $\Phi (1)$, $\Phi (2)$, and $\Phi (3)$ are the standard deviations of the critical source problem.

The FOM (a figure of merit) is calculated for one tally bin of each tally as a function of the sample numbers. It is defined as: $FOM=1/({\sigma }^{2}t)$ (2) where t is the computation time (min). FOM value was used to determine the Monte Carlo calculation efficiency. The larger FOM means the lesser the computer time t and a small variance s.

Variance reduction technologies

Reducing the variance σ is important for studying the Monte Carlo method. In JMCT, the main techniques include weight window, importance, mesh window, and source biasing. Biasing parameters may be produced based on adjoint calculations. At present, JMCT has a function of adjoint calculation. Other variance reduction techniques include geometry splitting and Russian roulette, energy splitting/roulette, implicit capture, forced collisions, source variable biasing, correlated sampling, exponential transformation, time cutoff, weight cutoff and energy cutoff, direction probability for point detector estimator, and DXTRAN etc.

Tallies

All tallies can express as: $I=\langle \varphi ,g\rangle ={\displaystyle {\int }_D\varphi (P)g(P)\text{d}P}$ (3) where $P=(r,E,\text{Ω},t)$, $\varphi (P)$ is a flux, $g(P)$ is a responding function.

$\varphi (P)$ satisfies the Boltzmann differential-integral equation as follows ^[20]: $\begin{array}{l}\frac{\text{1}}{\text{v}}\frac{\partial \varphi (r,E,\text{Ω},t)}{\partial t}+\text{Ω}\cdot \nabla \varphi (r,E,\text{Ω},t)+{\Sigma }_{\text{t}}\varphi \\ ={\displaystyle \iint {\Sigma }_{\text{s}}(r,E\text{'},\text{Ω}\text{'}\to E,\text{Ω})}\varphi (r,E\text{'},\text{Ω}\text{'},t)d{E}^{\prime }d\text{Ω}\text{'}\\ \text{+ }\frac{\chi (r,E)}{4\pi k_{\text{eff}}}{\displaystyle \iint \nu {\Sigma }_{\text{f}}(r,E)}\varphi (r,E\text{'},\text{Ω}\text{'},t)\text{d}{E}^{\prime }d\text{Ω}\text{'}+S(r,E,\text{Ω},t).\end{array}$ (4)

In the JMCT, the standard tallies include k-effective, current, surface flux, volume flux, point flux, deposition energy, detector response, and various reactive rates. In addition, the JMCT also provides a standard port of tally subroutine for user.

Domain decomposition and parallelization [21]

Currently, most Monte Carlo particle transport programs do not have a DD function which includes the MCNP program. The Mercury Monte Carlo particle transport program was the first program with a DD function ^[22]. If the case happens when one zone is divided into two or more zones (see Fig. 6(a)), there are some differences between a DD case and the no DD case. For JMCT another DD algorithm has been developed, in which the nature of the geometry interface is determined by random sampling (see Fig. 6(b)). Due to the topology structure maintaining, the same results are obtained in the DD case and no DD case. Of cause, it is necessary for designing a clear chain of random numbers (see Fig. 7).

Fig. 6Full-screen View

(Color online) Result comparison of DD and without DD case. (a) DD schematic of mercury; (b) DD schematic of JMCT.

Fig. 7Full-screen View

(Color online) Pin flux comparison of a reactor core in DD and no DD. (a) Pin flux without DD; (b) Pin flux with 2×1×1 DD; (c) Pin flux with 2×2×1 DD.

JMCT supports two-level parallelization of MPI and OpenMP. Figure 8(a) shows the communication between different domains in the DD case. Figure 8(b) shows domain replication, domain decomposition, and parallelization. Dr. Gang Li developed several algorithms to solve data exchange, dynamic load balance, and variance lower questions ^[21]. Figure 8(c) shows a comparison of the parallel efficiency in the DD and no DD case. A high parallel efficiency was still obtained even in the DD case.

Fig. 8Full-screen View

Domain communication and comparison of parallel efficiency in DD and no DD. (a) Domain communication in 2×2 case; (b) Domain replication and decomposition; (c) Parallel efficiency for 2×1DD, 2×2DD and no DD.

Uniform tally density (UTD)

Monte Carlo methods are widely used for criticality calculations. For this purpose, the power iteration method is one of the most important technique for the Monte Carlo program. To obtain an accurate multiplication factor and global tallies, inactive iteration cycles must be performed until the source distribution converges. Generally, tallying should only be invoked in the subsequent active cycles after completion of the inactive cycles. If the directed simulation uses, the errors in some margin assemblies may exceed the convergence rule of 95%. Therefore, the MC21 Monte Carlo program proposes a uniform fission site (UFS) algorithm for computing the expected number of fission secondary neutrons. UFS algorithm significantly improves the convergence of margin assemblies ^[23]. For JMCT, an improved algorithm based on the UFS algorithm has been proposed. This algorithm is called the uniform tally density (UTD) algorithm, and the modified formula is as follows: $m_{\text{JMCT}}^{\text{std}}=w_{\text{in}}\frac{\nu {\Sigma }_{\text{f}}}{{\widehat{k}}_{\text{eff}}{\Sigma }_t}\frac{v_k}{t_k}$ (5) where m is the multiplication factor estimated in the last cycle, v_k is a fraction of V occupied by zone k, in which the collision occurs, V is the volume of the problem domain that comprises all fissionable materials, s_k is a fraction of the fission source contained in zone k. For same problem, a large FOM value is obtained based on UTD ^[24].

Doppler broadening on-the-fly (OTF)

The on-the-fly (OTF) Doppler method was used to calculate the temperature-dependent cross-section. The OTF method can be used to calculate the cross-section of any nuclide at any temperature in the range of 300-3000 K based on a cross-section library of 300 K. For the OTF method, first, a series of temperature-dependent cross-section is produced by NJOY ^[11]. Second, a uniform energy grid was evaluated using temperature-dependent cross-section. Third, a polynomial was used to fit the temperature-dependent cross section of each energy grid. The coefficient of the polynomial was obtained using a single-value decomposition algorithm. The basic formula is as follows: $\begin{array}{c}{\sigma }_x(v,T)={\displaystyle \sum _{i=1}^n\frac{a_i}{T^{i/2}}+}{\displaystyle \sum _{i=1}^n{b}_i{T}^{i/2}+C},T\\ =300~3000\text{K},\Delta T=25\text{K};x\text{=}\text{t},\text{a},\text{f}\end{array}$ (6) where t, a, f is the total, absorption, and fission cross sections. Finally, the coefficients of the polynomial in all energy grids and the energy grids themselves were written in a text file for interpolation ^[25].

Fast criticality search of boron concentration (FCSBC)

The algorithms are based on the neutron balance equation. When considering significant reactions, such as (n,a), (n,f), (n,γ), and (n,2n), the following equation is obtained: $\begin{array}{l}{\displaystyle {\int }_v\left({\rho }_{\text{b}}{\sigma }_{\text{boron}}+{\Sigma }_{n,\gamma }+{\Sigma }_{\text{f}}+\text{Ω}\cdot \nabla -{\Sigma }_{n,2n}\right)}\varphi \text{(}\text{r})\text{d}\text{r}\\ =\frac{1}{k_{\text{eff}}}{\displaystyle {\int }_{v}F\varphi \text{(}\text{r})\text{d}\text{r}}\end{array}$ (7) where ${\sigma }_{\text{boron}}$ is a microscopic cross-section of soluble boron, ∑ are the related macroscopic cross-sections, F present the neutron produced by fission. Assuming that the ratio of soluble boron density to water density is the same in the entire reactor, the boron concentration can be iterated together with $k_{\text{eff}}$ during the cycles. The new algorithm is only one cycle iteration for convergence and the traditional method usually needs four to five time iterations ^[26].

Coupled MC and SN

It is well known that CADIS (or FW-CADIS) is an ideal method for solving the deep penetration problem. The importance values of Monte Carlo transport are produced by S_N code. However, it needs two modeling for S_N and MC calculation, respectively. It is good luck that JMCT (MC) and JSNT (S_N) use same modeling, it makes the CADIS method become easy and highly efficient. The mesh weight windows, importance values, and source biasing coefficients were produced by solving the adjoint equation based on the discrete ordinate S_N program JSNT ^[12]. Figure 9 shows the coupled flow of MC and S_N. Good efficiency has been achieved for some deep-penetration shielding problems ^{[27, 28]}.

Fig. 9Full-screen View

Flow of the coupled MC and S_N

Multi-physics coupled calculation

Tightly coupled neutron transport (JMCT), depletion (JBURN), and thermal hydraulics (JTH) have been developed for numerical reactor simulations. Recently, the JMCT has added the analysis of fuel properties and makes the simulation closer to the true. Figure 10 shows the iteration flow.

Fig. 10Full-screen View

Iterate flow of neutronics, thermal hydraulics and fuel property. (a) Data exchange among neutronics, thermal-hydraulic, and fuel properties; (b) Flow diagram of the multiphysics calculation.

Collision mechanism based on material [29]

For most Monte Carlo multigroup neutron transport calculations, the collision mechanism is based on nuclides. A new collision mechanism has been developed based on the material. This mechanism is suitable for the case of a large amount of fission production. The simulation can save a significant amount of memory and computational time. The basic principle is that suppose the material k being consists of m(k) types of nuclides. The combined cross sections of material k are defined as follows: ${\sigma }_x^g(k)={\displaystyle \sum _{i=1}^{m(k)}{n}_i{\sigma }_{x,i}^g}{,}_{}x=\text{t},\text{a},\text{f}$ (8) where n_i (i = 1,2, ..., m(k)) is the percentage of each nuclide, and g, t, a, f are the energy group, total, absorption, and fission cross-sections, respectively. The scattering cross section is defined as follows: ${\sigma }_s^{g\text{'}\to g}(k)={\displaystyle \sum _{i=1}^{m(k)}{n}_i{\sigma }_{s,i}^{g\text{'}\to g}}$ (9)

If ${\sigma }_{\text{f}}^g\ne 0$, then the number of fission neutrons is defined as follows: ${\nu }_g(k)={(\nu \sigma )}_{\text{f}}^g(k)/{\sigma }_{\text{f}}^g(k)$ (10)

The fission spectrum is defined as follows: ${\chi }_g(k)={\displaystyle \sum _{i=1}^{m(k)}{n}_i{(\nu \sigma )}_{\text{f},i}^g{\chi }_{g,i}}/{\displaystyle \sum _{i=1}^{m(k)}{n}_i{(\nu \sigma )}_{\text{f},i}^g}$ (11) where ${\chi }_{g,i}$ is the fission spectrum of i fission nuclide; g = 1,2,…,G, where g=1 corresponds to the highest energy group and g = G corresponds to the lowest energy group.

The scattering transfer matrix is defined as follows: $\left[\begin{array}{cccc}{\sigma }_{\text{s}}^{1\to 1}& {\sigma }_{\text{s}}^{1\to 2}& {\sigma }_{\text{s}}^{1\to 3}& \begin{array}{cc}\cdots & {\sigma }_{\text{s}}^{1\to G}\end{array}\\ {\sigma }_{\text{s}}^{2\to 1}& {\sigma }_{\text{s}}^{2\to 2}& {\sigma }_{\text{s}}^{2\to 2}& \begin{array}{cc}\cdots & {\sigma }_{\text{s}}^{2\to G}\end{array}\\ {\sigma }_{\text{s}}^{3\to 1}& {\sigma }_{\text{s}}^{3\to 2}& {\sigma }_{\text{s}}^{3\to 3}& \begin{array}{cc}\cdots & {\sigma }_{\text{s}}^{3\to G}\end{array}\\ \begin{array}{c}\cdots \\ 0\end{array}& \begin{array}{c}\cdots \\ \cdots \end{array}& \begin{array}{c}\cdots \\ {\sigma }_{\text{s}}^{G\to G-\text{nup}}\end{array}& \begin{array}{c}\begin{array}{cc}\cdots & \cdots \end{array}\\ \begin{array}{cc}\cdots & {\sigma }_{\text{s}}^{G\to G}\end{array}\end{array}\end{array}\right]$ (12)

If upper scattering is not considered, that is, nup = 0, then the scattering transfer matrix becomes a triangular matrix.

Treatment of multigroup angular distribution

The multigroup scattering cross section can be expressed as follows: ${\sigma }_{\text{s}}^{g\text{'}\to g}(x,\mu )={\sigma }_{\text{s}}^{g\text{'}\to g}(x)f(\mu )$ (13) where $\mu =\text{Ω}\text{'}\cdot \text{Ω}$; $\text{Ω}\text{'}$, and $\text{Ω}$ are the incident and emission directions, respectively, and $f(\mu )$ is an angle distribution function that satisfies normalization, that is, $f(\mu )\ge 0{,}_{}{\displaystyle {\int }_{-1}^{1}f(\mu )d\mu }=1$.

In general, $f(\mu )$ is expanded in a Legendre series as following $f(\mu )={\displaystyle \sum _{l=0}^{\infty }\frac{2l+1}{2}{f}_l{P}_l(\mu )}$ (14) where Legendre coefficients are: $f_l={\sigma }_{s,l}^{g\text{'}\to g}/{\sigma }_{s,0}^{g\text{'}\to g}$ (15)

The approximate distribution is obtained by making L-order cutoff: $f_L(\mu )={\displaystyle \sum _{l=0}^L\frac{2l+1}{2}{f}_l{P}_l(\mu )}$ (16)

Due to this cutoff, $f_L(\mu )$ may appear negative in [-1,1]. To avoid this case, the generalized Gaussian quadrature is used ^[30] and an approximate discrete distribution is obtained as follows: $f^{*}(\mu )={\displaystyle \sum _{i=1}^n{p}_i\delta (\mu -{\mu }_i)}$ (17) where $n=\frac{L+1}{2}$, $p_i={\left[{\displaystyle \sum _{k=1}^{n-1}{Q}_k^2({\mu }_i)}/N_k\right]}^{-1}$ is the probability of the sample ${\mu }_i$ (${\displaystyle {\sum }_{i=1}^n{p}_i}=1$). ${\left\{Q_i(\mu )\right\}}_{i=1}^n$ is an orthogonal polynomial system with respect to $f_L(\mu )$, ${\left\{{\mu }_i\right\}}_{i=1}^n$ is the root of $Q_n(\mu )$ and satisfies: ${\displaystyle {\int }_{-1}^1{Q}_i(\mu )}{Q}_j(\mu )f_L(\mu )\text{d}\mu ={\delta }_{ij}{N}_i$ (18)

where ${\delta }_{ij}$ is the Kronecker delta function and $N_i$ is a normalization constant ^[29]. $f^{*}(\mu )$ replaces $f_L(\mu )$ and m is produced through sampling.

Random geometry modeling

The JMCT is used to simulate a high-temperature reactor (HTR). The core of the HTR consists of many fuel pebbles, with each fuel pebble consisting of many individual pellets ^[31]. According to the traditional method, these pellets are uniformly mixed into a single material. Today, precision physical research is needed to study the effect of pellet numbers and different distributions on pebble criticality. Therefore, the geometry must describe according to the real distribution of pellet. So, a sample method has been developed based on the real distribution of pellet. Figure 11 shows the uniform pellet distribution in the pebble.

Fig. 11Full-screen View

Pellet distribution in pebble based on the sampling method

HZP condition

As a marked achievement of the CASL plan ^[34], VERA core physics benchmark progression problem specifications were published on March 29, 2013 ^[5]. Figure 17 shows the VERA modeling by the JLAMT. Table 4 compares the k_eff and control rod worth between the JMCT and KENO-VI ^[10]. Figure 18 shows a comparison of the radial integral power distributions of the JMCT and KENO-VI. The maximum difference was 0.85%. Good agreement was achieved between the JMCT and the KENO-VI.

Fig. 17Full-screen View

(Color online) VERA geometry by JMCT

Fig. 18Full-screen View

Comparison of assembly power between JMCT and KENO-VI. (a) Radial assembly power distribution; (b) Axial relative power distribution.

MAX=0.85%RMS=0.36%

HFP condition

Figure 19 shows a comparison of the radial integral assembly power distribution between the JMCT and MC21. The maximum difference was 2.27%.

Fig. 19Full-screen View

Comparison of the radial integral assembly power distribution between JMCT and MC21

Reactor core

For CAP1400, the initial physical parameters of the first core were calculated in a similar fashion to that of AP1000. The JMCT simulation results were selected as the reference values. Table 6 lists the deviations of NECP-X and SCAP relative to JMCT ^[37].

Shielding

The radiation shielding of the CAP1400 includes eight operating cases (i.e., eight models). Here, only the simulation results for a typical case are presented. First, a complex source subroutine is designed according to the distribution of position, energy, and direction. The fixed source and mesh are tally adopted, where the importance values are produced by the S_N program JSNT. 26000 million neutron histories were simulated using 260000 processors. The total CPU time was 17.7 minutes. The calculation was completed using a Tianhe-II computer. Figure 22 shows the neutron dose distribution and secondary photon dose distribution of the full factory.

Fig. 22Full-screen View

(Color online) Model and simulation results of CAP1400 full factory. (a) Shield model of CAP1400; (b) Neutron dose distribution; (c) Secondary photon dose distribution.