A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0 Cover Article Enhanced Publication

A high fidelity general purpose 3-D Monte Carlo particle transport program JMCT3.0 1 2 first-author https://orcid.org/0000-0002-5786-3603 1 2 corresp *li_gang@iapcm.ac.cn *li_gang@iapcm.ac.cn li_gang@iapcm.ac.cn 1 2 2 2 2 https://orcid.org/0000-0002-4118-4075 2 2 2 1 2 1 1 2 https://orcid.org/0000-0001-7585-5167 2 Institute of Applied Physics and Computational Mathematics (IAPCM), Beijing 100094, China Institute of Applied Physics and Computational Mathematics (IAPCM), Beijing 100094, China CAEP Software Centre for High Performance Numerical Simulation (CAEP-SCNS), Beijing 100088, China CAEP Software Centre for High Performance Numerical Simulation (CAEP-SCNS), Beijing 100088, China JMCT is a large-scale, high-fidelity, three-dimensional general neutron-photon-electron-proton transport Monte Carlo software system. It was developed based on the combinatorial geometry parallel infrastructure JCOGIN and the adaptive structured mesh infrastructure JASMIN. JMCT is equipped with CAD modeling and visualizes the image output. It supports the geometry of the body and the structured/unstructured mesh. JMCT has most functions, variance reduction techniques, and tallies of the traditional Monte Carlo particle transport codes. Two energy models, multi-group and continuous, are provided. In recent years, some new functions and algorithms have been developed, such as Doppler broadening on-the-fly (OTF), uniform tally density (UTD), consistent adjoint driven importance sampling (CADIS), fast criticality search of boron concentration (FCSBC), domain decomposition (DD), adaptive control rod moving (ACRM), and random geometry (RG) etc. The JMCT is also coupled with the discrete ordinate S N code JSNT to generate source-biasing factors and weight-window parameters. At present, the number of geometric bodies, materials, tallies, depletion zones, and parallel processors are sufficiently large to simulate extremely complicated device problems. JMCT can be used to simulate reactor physics, criticality safety analysis, radiation shielding, detector response, nuclear well logging, and dosimetry calculations. In particular, JMCT can be coupled with depletion and thermal-hydraulics for the simulation of reactor nuclear-hot feedback effects. This paper describes the progress in advanced modeling, high-performance numerical simulation of particle transport, multiphysics coupled calculations, and large-scale parallel computing. Advanced modeling High performance numerical simulation Multi-physics coupled calculation Large-scale parallel computing JMCT Introduction 1 Today, numerical simulations allow us to obtain as much detailed information as possible for any nuclear system based on the Monte Carlo (MC) method and software. The physical quantities can easily obtain both measurable and immeasurable values, especially for some extreme conditions inaccessible by experiments. Monte Carlo methods and software have been widely used in nuclear science engineering, statistical physics, biomedicine, quantum mechanics, molecular dynamics, petroleum geophysical exploration, finance, information, operational research, and polymer chemistry. JMCT (Joint Monte Carlo Transport) has been developed by the IAPCM and CAEP-SCNS. It is designed to simulate neutrons, photons, electrons, proton, light radiation, atmosphere and molecule transport. This passes the validation and verification of a large number of benchmarks and experiments. The simulation problems cover the standard, critical fission, and adjoint sources. The data libraries use the continuous point-wise cross-section of the ACE format and multigroup cross-section of the ANISN format. Abundant variance reduction techniques have been incorporated into the JMCT. The JMCT software has been developed based on the parallel infrastructure JCOGIN [ 1 1 ] and JASMIN [ 2 2 ] , where JCOGIN is a parallel combinatorial geometry infrastructure and JASMIN is an adaptive structured mesh infrastructure. The input uses CAD modeling, and the calculated results are outputted as a visual format or text file. Some problems in which the memory exceeds the limit of a single core or node can be simulated using the DD algorithm. At present, the number of geometric bodies, materials, tallies, depletion zones, sample numbers, period of random numbers, and parallel processors are large enough and expandable. In particular, the JMCT is coupled depletion and thermal hydraulic, where the depletion supports the TTA and CRAM methods, and thermal hydraulic supports the sub-channel method [ 3 3 ] . In this paper, first, the basic features and key algorithms are introduced. Then the numerical results of the BEAVRS [ 4 4 ] , VERA [ 5 5 ] , AP1000 [ 6 6 ] and CAP1400 [ 7 7 ] reactor models are presented. For the BEAVRS, VERA, and AP1000 models, the numerical results of the JMCT are in good agreement with the results of MC21 [ 8 8 ] , OpenMC [ 9 9 ] , KENO-VI [ 10 10 ] and experiments. Finally, the initial physical parameters of the CAP1400 first core were predicted. The high accuracy of the JMCT was demonstrated by comparison with other well-known Monte Carlo programs. JMCT and methods 1 JMCT system 2 JMCT is a program of the JPTS (Joint Particle Transport System) package. The JPTS package includes three floors, where the top floor contains the data library NuDa produced by NJOY [ 11 11 ] . The middle floor contains four application programs: JSNT (S N ) [ 12 12 ] , JMCT [ 3 3 ] (MC), JBURN [ 13 13 ] (deletion), and JTH (thermal-hydraulic). The lower floor comprises the two support frameworks of JASMIN [ 2 2 ] and JCOGIN [ 1 1 ] . The unified CAD preprocessor JLAMT [ 14 14 ] and view postprocessor TeraVAP [ 15 15 ] are equipped for the JPTS (see Fig. 1 Fig. 1 ). Figure 2 Figure 2 shows the structure flow of the JMCT, where the JMCT is coupled with the JSNT for radiation shielding and the JMCT is coupled with the JBURN and JTH for the reactor core. JMCT modeling and geometry 2 The JMCT input is performed with the help of preprocess software JLAMT, which provides CAD modeling. Figure 3 Figure 3 shows the input interface. The input parameters include geometry, material, density, temperature, importance, tally types, source information, sample numbers, a wide variety of variance reduction techniques, energy cutoff, weight cutoff, and time cutoff etc. The material in the geometric zone consists of multiple nuclei, including isotopes. Specifically, JMCT supports repetitive structures of the same geometry zone with different materials (MCNP [ 16 16 ] only supports repetitive structures of the same geometry zone with the same material). This function is available for the tightly coupled of neutron transport and depletion. The JMCT tally assorts a global tally and user-specified tally. By default, a global tally is performed for all zones according to the track length estimator. In addition, a mesh tally was added to the structured mesh. The JMCT was developed on the basis of the JCOGIN infrastructure. JCOGIN supports the real combination geometry body and it integrates the common Monte Carlo particle transport methods, such as geometry descriptions, Boolean operations, track length calculations, random number generators, sample of arbitrary distributions, calculation of volumes and areas, domain decomposition, parallel computations, and load balance technologies etc. At present, the JCOGIN infrastructure has supported the development of several Monte Carlo programs, and the JMCT is one of these programs. JMCT uses real combinatorial geometry, support geometry bodies, and structured/unstructured meshes. The basic geometric bodies include cube, cylinder, sphere, cone, ellipsoid, and torus etc. Some special bodies, such as reactor components, 1/n pillars, 1/n spheres, and special fuel assemblies, can be customized (see Fig. 4 Fig. 4 ). JMCT output 2 The JMCT supports visual output and text file output. The postprocessor TeraVAP supports image output in TBs scale. Figure 5 Figure 5 shows some examples. JMCT features 2 Nuclear data 3 The data were from the ENDF/B-VII [ 17 17 ] and CENDL-3.2 [ 18 18 ] evaluation libraries. The cross-section parameters were generated by using NJOY code [ 11 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 5 cross-section library is mainly used for radiation shielding, which its energy structure is same as the Bugle library [ 19 19 ] . Another 172-group neutrons, 40-group upper scatter neutrons, and 30-group photon P 5 library is mainly used for the reactor core. Source specification 3 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 3 Monte Carlo error e satisfies the center limit theorem as following <math display="block" id="M1"><mrow><munder><mrow><mi>lim</mi></mrow><mrow><mi>N</mi><mo>→</mo><mi>∞</mi></mrow></munder><mi>P</mi><mrow><mo>{</mo><mrow><mrow><mo>|</mo><mrow><mfrac><mrow><mover accent="true"><mi>x</mi><mo>¯</mo></mover><mo>−</mo><mi>a</mi></mrow><mover accent="true"><mi>x</mi><mo>¯</mo></mover></mfrac></mrow><mo>|</mo></mrow><mo><</mo><mi>ε</mi></mrow><mo>}</mo></mrow><mo>=</mo><mi>Φ</mi><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo></mrow></math> where <math id="M2"><mrow><mover accent="true"><mi>x</mi><mo>¯</mo></mover><mo>=</mo><mfrac><mn>1</mn><mi>N</mi></mfrac><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>N</mi></msubsup><mrow><msub><mi>x</mi><mi>i</mi></msub></mrow></mstyle></mrow></math> is the mean, N is the sample number (or histories), <math id="M3"><mrow><mi>a</mi><mo>=</mo><mi>E</mi><mo stretchy="false">(</mo><msub><mi>x</mi><mi>i</mi></msub><mo stretchy="false">)</mo><mo>,</mo><mi>i</mi><mo>=</mo><mn>1</mn><mo>,</mo><mn>2</mn><mo>,</mo><mo>⋯</mo></mrow></math> is a mathematical expectation of various random x i , <math id="M4"><mrow><mi>ε</mi><mo>=</mo><mfrac><mrow><mi>λ</mi><mi>σ</mi></mrow><mrow><msqrt><mi>N</mi></msqrt><mover accent="true"><mi>x</mi><mo>¯</mo></mover></mrow></mfrac></mrow></math> is an error, s σ is a variance and <math id="M5"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mi>λ</mi><mo stretchy="false">)</mo><mo>=</mo><mfrac><mn>1</mn><mrow><msqrt><mrow><mn>2</mn><mi>π</mi></mrow></msqrt></mrow></mfrac><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mi>λ</mi></mrow><mi>λ</mi></msubsup><mrow><msup><mi>e</mi><mrow><mo>−</mo><msup><mi>t</mi><mn>2</mn></msup><mo>/</mo><mn>2</mn></mrow></msup><mi>d</mi><mi>t</mi></mrow></mrow></mstyle></mrow></math> is the degree of confidence. <math id="M6"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo><mo>=</mo><mn>0.6827</mn></mrow></math> , <math id="M7"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>2</mn><mo stretchy="false">)</mo><mo>=</mo><mn>0.9545</mn></mrow></math> , <math id="M8"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>3</mn><mo stretchy="false">)</mo><mo>=</mo><mn>0.9973</mn></mrow></math> . In general, the error in <math id="M9"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> is the standard deviation of the fixed-source problem. The errors in <math id="M10"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>1</mn><mo stretchy="false">)</mo></mrow></math> , <math id="M11"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>2</mn><mo stretchy="false">)</mo></mrow></math> , and <math id="M12"><mrow><mi>Φ</mi><mo stretchy="false">(</mo><mn>3</mn><mo stretchy="false">)</mo></mrow></math> 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: <math display="block" id="M13"><mrow><mi>F</mi><mi>O</mi><mi>M</mi><mo>=</mo><mn>1</mn><mo>/</mo><mo stretchy="false">(</mo><msup><mi>σ</mi><mn>2</mn></msup><mi>t</mi><mo stretchy="false">)</mo></mrow></math> 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 3 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 3 All tallies can express as: <math display="block" id="M14"><mrow><mi>I</mi><mo>=</mo><mrow><mo>〈</mo><mrow><mi>ϕ</mi><mo>,</mo><mi>g</mi></mrow><mo>〉</mo></mrow><mo>=</mo><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>D</mi></msub><mrow><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">)</mo><mi>g</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle></mrow></mrow></mstyle></mrow></math> where <math id="M15"><mrow><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo>=</mo><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo>,</mo><mstyle mathvariant="bold" mathsize="normal"><mtext>Ω</mtext></mstyle><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo></mrow></math> , <math id="M16"><mrow><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">)</mo></mrow></math> is a flux, <math id="M17"><mrow><mi>g</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">)</mo></mrow></math> is a responding function. <math id="M18"><mrow><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>P</mi></mstyle><mo stretchy="false">)</mo></mrow></math> satisfies the Boltzmann differential-integral equation as follows [ 20 20 ] : <math display="block" id="M19"><mtable columnalign="left"><mtr><mtd><mfrac><mtext>1</mtext><mtext>v</mtext></mfrac><mfrac><mrow><mo>∂</mo><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold-italic" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo></mrow><mrow><mo>∂</mo><mi>t</mi></mrow></mfrac><mo>+</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>⋅</mo><mo>∇</mo><mi>ϕ</mi><mo stretchy="false">(</mo><mi>r</mi><mo>,</mo><mi>E</mi><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mo>+</mo><msub><mi>Σ</mi><mtext>t</mtext></msub><mi>ϕ</mi></mtd></mtr><mtr><mtd><mo>=</mo><mstyle displaystyle="true"><mrow><mo>∬</mo><mrow><msub><mi>Σ</mi><mtext>s</mtext></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo>'</mo><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>'</mo><mo>→</mo><mi>E</mi><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo stretchy="false">)</mo></mrow></mrow></mstyle><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo>'</mo><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>'</mo><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mi>d</mi><msup><mi>E</mi><mo>′</mo></msup><mi>d</mi><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>'</mo></mtd></mtr><mtr><mtd><mtext>+ </mtext><mfrac><mrow><mi>χ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo></mrow><mrow><mn>4</mn><mi>π</mi><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><mo>∬</mo><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo stretchy="false">)</mo></mrow></mrow></mstyle><mi>ϕ</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo>'</mo><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>'</mo><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mtext>d</mtext><msup><mi>E</mi><mo>′</mo></msup><mi>d</mi><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>'</mo><mo>+</mo><mi>S</mi><mo stretchy="false">(</mo><mstyle mathvariant="bold" mathsize="normal"><mi>r</mi></mstyle><mo>,</mo><mi>E</mi><mo>,</mo><mstyle mathvariant="bold"><mtext>Ω</mtext></mstyle><mo>,</mo><mi>t</mi><mo stretchy="false">)</mo><mo>.</mo></mtd></mtr></mtable></math> 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. Algorithms 2 In recent years, a series of special algorithms for different problems have been developed for JMCT. Here, only a few key algorithms are introduced. Domain decomposition and parallelization [ 21 21 ] 3 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 22 ] . If the case happens when one zone is divided into two or more zones (see Fig. 6(a) 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) 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. 7 ). JMCT supports two-level parallelization of MPI and OpenMP. Figure 8(a) Figure 8(a) shows the communication between different domains in the DD case. Figure 8(b) 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 21 ] . Figure 8(c) 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. Uniform tally density (UTD) 3 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 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: <math display="block" id="M20"><mrow><msubsup><mi>m</mi><mrow><mtext>JMCT</mtext></mrow><mrow><mtext>std</mtext></mrow></msubsup><mo>=</mo><msub><mi>w</mi><mrow><mtext>in</mtext></mrow></msub><mfrac><mrow><mi>ν</mi><msub><mi>Σ</mi><mtext>f</mtext></msub></mrow><mrow><msub><mover accent="true"><mi>k</mi><mo>^</mo></mover><mrow><mtext>eff</mtext></mrow></msub><msub><mi>Σ</mi><mi>t</mi></msub></mrow></mfrac><mfrac><mrow><msub><mi>v</mi><mi>k</mi></msub></mrow><mrow><msub><mi>t</mi><mi>k</mi></msub></mrow></mfrac></mrow></math> 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 24 ] . Doppler broadening on-the-fly (OTF) 3 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 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: <math display="block" id="M21"><mtable><mtr><mtd><msub><mi>σ</mi><mi>x</mi></msub><mo stretchy="false">(</mo><mi>v</mi><mo>,</mo><mi>T</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><mfrac><mrow><msub><mi>a</mi><mi>i</mi></msub></mrow><mrow><msup><mi>T</mi><mrow><mi>i</mi><mo>/</mo><mn>2</mn></mrow></msup></mrow></mfrac><mo>+</mo></mrow></mstyle><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><msub><mi>b</mi><mi>i</mi></msub><msup><mi>T</mi><mrow><mi>i</mi><mo>/</mo><mn>2</mn></mrow></msup><mo>+</mo><mi>C</mi></mrow></mstyle><mo>,</mo><mi>T</mi></mtd></mtr><mtr><mtd><mo>=</mo><mn>300</mn><mo>~</mo><mn>3000</mn><mtext> </mtext><mtext>K</mtext><mo>,</mo><mi>Δ</mi><mi>T</mi><mo>=</mo><mn>25</mn><mtext> </mtext><mtext>K</mtext><mo>;</mo><mi>x</mi><mtext>=</mtext><mtext>t</mtext><mo>,</mo><mi> </mi><mtext>a</mtext><mo>,</mo><mi> </mi><mtext>f</mtext></mtd></mtr></mtable></math> 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 25 ] . Fast criticality search of boron concentration (FCSBC) 3 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: <math display="block" id="M22"><mtable columnalign="left"><mtr><mtd><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>v</mi></msub><mrow><mrow><mo>(</mo><mrow><msub><mi>ρ</mi><mtext>b</mtext></msub><msub><mi>σ</mi><mrow><mtext>boron</mtext></mrow></msub><mo>+</mo><msub><mi>Σ</mi><mrow><mi>n</mi><mo>,</mo><mi>γ</mi></mrow></msub><mo>+</mo><msub><mi>Σ</mi><mtext>f</mtext></msub><mo>+</mo><mtext>Ω</mtext><mo>⋅</mo><mo>∇</mo><mo>−</mo><msub><mi>Σ</mi><mrow><mi>n</mi><mo>,</mo><mn>2</mn><mi>n</mi></mrow></msub></mrow><mo>)</mo></mrow></mrow></mrow></mstyle><mi>ϕ</mi><mtext>(</mtext><mstyle mathvariant="bold"><mtext>r</mtext></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold"><mtext>r</mtext></mstyle></mtd></mtr><mtr><mtd><mtext>  </mtext><mo>=</mo><mfrac><mn>1</mn><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><mrow><msub><mo>∫</mo><mi>v</mi></msub><mrow><mi>F</mi><mi>ϕ</mi><mtext>(</mtext><mstyle mathvariant="bold"><mtext>r</mtext></mstyle><mo stretchy="false">)</mo><mtext>d</mtext><mstyle mathvariant="bold"><mtext>r</mtext></mstyle></mrow></mrow></mstyle></mtd></mtr></mtable></math> where <math id="M23"><mrow><msub><mi>σ</mi><mrow><mtext>boron</mtext></mrow></msub></mrow></math> 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 <math id="M24"><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></math> 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 26 ] . Coupled MC and S N 3 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 12 ] . Figure 9 Figure 9 shows the coupled flow of MC and S N . Good efficiency has been achieved for some deep-penetration shielding problems [ 27 27 , 28 28 ] . Multi-physics coupled calculation 3 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 Figure 10 shows the iteration flow. Collision mechanism based on material [ 29 29 ] 3 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: <math display="block" id="M25"><mrow><msubsup><mi>σ</mi><mi>x</mi><mi>g</mi></msubsup><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>m</mi><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></munderover><mrow><msub><mi>n</mi><mi>i</mi></msub><msubsup><mi>σ</mi><mrow><mi>x</mi><mo>,</mo><mi>i</mi></mrow><mi>g</mi></msubsup></mrow></mstyle><msub><mo>,</mo><mrow/></msub><mi>x</mi><mo>=</mo><mtext>t</mtext><mo>,</mo><mtext>a</mtext><mo>,</mo><mtext>f</mtext></mrow></math> 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: <math display="block" id="M26"><mrow><msubsup><mi>σ</mi><mi>s</mi><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>m</mi><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></munderover><mrow><msub><mi>n</mi><mi>i</mi></msub><msubsup><mi>σ</mi><mrow><mi>s</mi><mo>,</mo><mi>i</mi></mrow><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup></mrow></mstyle></mrow></math> If <math id="M27"><mrow><msubsup><mi>σ</mi><mtext>f</mtext><mi>g</mi></msubsup><mo>≠</mo><mn>0</mn></mrow></math> , then the number of fission neutrons is defined as follows: <math display="block" id="M28"><mrow><msub><mi>ν</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo><mo>=</mo><msubsup><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>σ</mi><mo stretchy="false">)</mo></mrow><mtext>f</mtext><mi>g</mi></msubsup><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo><mo>/</mo><msubsup><mi>σ</mi><mtext>f</mtext><mi>g</mi></msubsup><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></math> The fission spectrum is defined as follows: <math display="block" id="M29"><mrow><msub><mi>χ</mi><mi>g</mi></msub><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>m</mi><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></munderover><mrow><msub><mi>n</mi><mi>i</mi></msub><msubsup><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mo>,</mo><mi>i</mi></mrow><mi>g</mi></msubsup><msub><mi>χ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></mstyle><mo>/</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>m</mi><mo stretchy="false">(</mo><mi>k</mi><mo stretchy="false">)</mo></mrow></munderover><mrow><msub><mi>n</mi><mi>i</mi></msub><msubsup><mrow><mo stretchy="false">(</mo><mi>ν</mi><mi>σ</mi><mo stretchy="false">)</mo></mrow><mrow><mtext>f</mtext><mo>,</mo><mi>i</mi></mrow><mi>g</mi></msubsup></mrow></mstyle></mrow></math> where <math id="M30"><mrow><msub><mi>χ</mi><mrow><mi>g</mi><mo>,</mo><mi>i</mi></mrow></msub></mrow></math> 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: <math display="block" id="M31"><mrow><mrow><mo>[</mo><mrow><mtable><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>1</mn><mo>→</mo><mn>1</mn></mrow></msubsup></mrow></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>1</mn><mo>→</mo><mn>2</mn></mrow></msubsup></mrow></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>1</mn><mo>→</mo><mn>3</mn></mrow></msubsup></mrow></mtd><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>1</mn><mo>→</mo><mi>G</mi></mrow></msubsup></mrow></mtd></mtr></mtable></mrow></mtd></mtr><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>2</mn><mo>→</mo><mn>1</mn></mrow></msubsup></mrow></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>2</mn><mo>→</mo><mn>2</mn></mrow></msubsup></mrow></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>2</mn><mo>→</mo><mn>2</mn></mrow></msubsup></mrow></mtd><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>2</mn><mo>→</mo><mi>G</mi></mrow></msubsup></mrow></mtd></mtr></mtable></mrow></mtd></mtr><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>3</mn><mo>→</mo><mn>1</mn></mrow></msubsup></mrow></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>3</mn><mo>→</mo><mn>2</mn></mrow></msubsup></mrow></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>3</mn><mo>→</mo><mn>3</mn></mrow></msubsup></mrow></mtd><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mn>3</mn><mo>→</mo><mi>G</mi></mrow></msubsup></mrow></mtd></mtr></mtable></mrow></mtd></mtr><mtr><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable></mrow></mtd><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd></mtr><mtr><mtd><mo>⋯</mo></mtd></mtr></mtable></mrow></mtd><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd></mtr><mtr><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mi>G</mi><mo>→</mo><mi>G</mi><mo>−</mo><mtext>nup</mtext></mrow></msubsup></mrow></mtd></mtr></mtable></mrow></mtd><mtd><mrow><mtable><mtr><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd><mtd><mo>⋯</mo></mtd></mtr></mtable></mrow></mtd></mtr><mtr><mtd><mrow><mtable><mtr><mtd><mo>⋯</mo></mtd><mtd><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mi>G</mi><mo>→</mo><mi>G</mi></mrow></msubsup></mrow></mtd></mtr></mtable></mrow></mtd></mtr></mtable></mrow></mtd></mtr></mtable></mrow><mo>]</mo></mrow></mrow></math> If upper scattering is not considered, that is, nup = 0, then the scattering transfer matrix becomes a triangular matrix. Treatment of multigroup angular distribution 3 The multigroup scattering cross section can be expressed as follows: <math display="block" id="M32"><mrow><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo>,</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><msubsup><mi>σ</mi><mtext>s</mtext><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup><mo stretchy="false">(</mo><mi>x</mi><mo stretchy="false">)</mo><mi>f</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> where <math id="M33"><mrow><mi>μ</mi><mo>=</mo><mstyle mathvariant="bold" mathsize="normal"><mtext>Ω</mtext></mstyle><mo>'</mo><mo>⋅</mo><mstyle mathvariant="bold" mathsize="normal"><mtext>Ω</mtext></mstyle></mrow></math> ; <math id="M34"><mrow><mstyle mathvariant="bold" mathsize="normal"><mtext>Ω</mtext></mstyle><mo>'</mo></mrow></math> , and <math id="M35"><mstyle mathvariant="bold" mathsize="normal"><mtext>Ω</mtext></mstyle></math> are the incident and emission directions, respectively, and <math id="M36"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> is an angle distribution function that satisfies normalization, that is, <math id="M37"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>≥</mo><mn>0</mn><msub><mo>,</mo><mrow/></msub><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mi>d</mi><mi>μ</mi></mrow></mrow></mstyle><mo>=</mo><mn>1</mn></mrow></math> . In general, <math id="M38"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> is expanded in a Legendre series as following <math display="block" id="M39"><mrow><mi>f</mi><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>l</mi><mo>=</mo><mn>0</mn></mrow><mi>∞</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msub><mi>f</mi><mi>l</mi></msub><msub><mi>P</mi><mi>l</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></math> where Legendre coefficients are: <math display="block" id="M40"><mrow><msub><mi>f</mi><mi>l</mi></msub><mo>=</mo><msubsup><mi>σ</mi><mrow><mi>s</mi><mo>,</mo><mi>l</mi></mrow><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup><mo>/</mo><msubsup><mi>σ</mi><mrow><mi>s</mi><mo>,</mo><mn>0</mn></mrow><mrow><mi>g</mi><mo>'</mo><mo>→</mo><mi>g</mi></mrow></msubsup></mrow></math> The approximate distribution is obtained by making L -order cutoff: <math display="block" id="M41"><mrow><msub><mi>f</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>l</mi><mo>=</mo><mn>0</mn></mrow><mi>L</mi></munderover><mrow><mfrac><mrow><mn>2</mn><mi>l</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac><msub><mi>f</mi><mi>l</mi></msub><msub><mi>P</mi><mi>l</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></mstyle></mrow></math> Due to this cutoff, <math id="M42"><mrow><msub><mi>f</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> may appear negative in [-1,1]. To avoid this case, the generalized Gaussian quadrature is used [ 30 30 ] and an approximate discrete distribution is obtained as follows: <math display="block" id="M43"><mrow><msup><mi>f</mi><mo>*</mo></msup><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></munderover><mrow><msub><mi>p</mi><mi>i</mi></msub><mi>δ</mi><mo stretchy="false">(</mo><mi>μ</mi><mo>−</mo><msub><mi>μ</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mstyle></mrow></math> where <math id="M44"><mrow><mi>n</mi><mo>=</mo><mfrac><mrow><mi>L</mi><mo>+</mo><mn>1</mn></mrow><mn>2</mn></mfrac></mrow></math> , <math id="M45"><mrow><msub><mi>p</mi><mi>i</mi></msub><mo>=</mo><msup><mrow><mrow><mo>[</mo><mrow><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><mi>k</mi><mo>=</mo><mn>1</mn></mrow><mrow><mi>n</mi><mo>−</mo><mn>1</mn></mrow></munderover><mrow><msubsup><mi>Q</mi><mi>k</mi><mn>2</mn></msubsup><mo stretchy="false">(</mo><msub><mi>μ</mi><mi>i</mi></msub><mo stretchy="false">)</mo></mrow></mstyle><mo>/</mo><msub><mi>N</mi><mi>k</mi></msub></mrow><mo>]</mo></mrow></mrow><mrow><mo>−</mo><mn>1</mn></mrow></msup></mrow></math> is the probability of the sample <math id="M46"><mrow><msub><mi>μ</mi><mi>i</mi></msub></mrow></math> ( <math id="M47"><mrow><mstyle displaystyle="true"><msubsup><mo>∑</mo><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup><mrow><msub><mi>p</mi><mi>i</mi></msub></mrow></mstyle><mo>=</mo><mn>1</mn></mrow></math> ). <math id="M48"><mrow><msubsup><mrow><mrow><mo>{</mo><mrow><msub><mi>Q</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow><mo>}</mo></mrow></mrow><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup></mrow></math> is an orthogonal polynomial system with respect to <math id="M49"><mrow><msub><mi>f</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> , <math id="M50"><mrow><msubsup><mrow><mrow><mo>{</mo><mrow><msub><mi>μ</mi><mi>i</mi></msub></mrow><mo>}</mo></mrow></mrow><mrow><mi>i</mi><mo>=</mo><mn>1</mn></mrow><mi>n</mi></msubsup></mrow></math> is the root of <math id="M51"><mrow><msub><mi>Q</mi><mi>n</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> and satisfies: <math display="block" id="M52"><mrow><mstyle displaystyle="true"><mrow><msubsup><mo>∫</mo><mrow><mo>−</mo><mn>1</mn></mrow><mn>1</mn></msubsup><mrow><msub><mi>Q</mi><mi>i</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></mrow></mstyle><msub><mi>Q</mi><mi>j</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><msub><mi>f</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo><mtext>d</mtext><mi>μ</mi><mo>=</mo><msub><mi>δ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub><msub><mi>N</mi><mi>i</mi></msub></mrow></math> where <math id="M53"><mrow><msub><mi>δ</mi><mrow><mi>i</mi><mi>j</mi></mrow></msub></mrow></math> is the Kronecker delta function and <math id="M54"><mrow><msub><mi>N</mi><mi>i</mi></msub></mrow></math> is a normalization constant [ 29 29 ] . <math id="M55"><mrow><msup><mi>f</mi><mo>*</mo></msup><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> replaces <math id="M56"><mrow><msub><mi>f</mi><mi>L</mi></msub><mo stretchy="false">(</mo><mi>μ</mi><mo stretchy="false">)</mo></mrow></math> and m is produced through sampling. Random geometry modeling 3 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 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 Figure 11 shows the uniform pellet distribution in the pebble. Numerical tests 1 BEAVRS model 2 The BEAVRS model is thought to be a significant challenge for computers and software. The model was released by the MIT Computational Reactor Physics Group on July 7, 2013 (www.crpg.mit.edu) [ 4 4 ] . It includes detailed specifications of the operating 4-loop Westinghouse PWR (3411 MW), two cycles of measured data, hot zero power (HZP) and hot full power (HFP) data, fuel loads by assembling as built, and three enriched fuels (1.6%, 2.4%, and 3.1%). Two cycles of measured data were used to validate the high-fidelity core analysis program. The simulation results of MC21 and OpenMC in the HZP condition were presented in PHYSOR2014 [ 32 32 ] . The basic parameters are as follows: (1) Fuel assembles: 193 (2) Axial planes: 398 (3) Pins/assemble:289 (17 × 17, 264 fuel pins, and 25 guide tubes) (4) Total tally regions: 193 × 17 × 17 × 1 × 398 = 22199246 (5) Total regions: 193 × 17 × 17 × 2 × 398 = 44398492 (6) Convergence standard: standard deviation of 95% pin powers ≤ 1% The core contained eight types of rods, as shown in Fig. 12(a) Fig. 12(a) , and nine types of assemblies, as shown in Fig. 12(b) Fig. 12(b) . Each assembly had six segment grids in the axial direction, as shown in Fig. 12(c) Fig. 12(c) . Because the model memory is too large, 2×´2×´2 DD are used, as shown in Fig. 12(d) Fig. 12(d) . The simulation tracked 4 million neutrons per cycle, 1000 cycles in total, and discarded 400 cycles. The simulation consumed 200 processors and 5.3 hours in the Tianhe-II computer. Figure 13 Figure 13 shows a comparison of the radial integral assembly power between the MC21 and JMCT. The maximal difference is 3.17% in the E2 assembly, and the minimum difference is 0 in the E6 assembly (MC21 appears in the same assemblies). Figure 14 Figure 14 shows a comparison of the detector tallies and measured data in the instrument tubes. The maximal deviation of the power assemble was: -14.77% in the B13 assembly (MC21 was: -17.05%). The minimum deviation of the power assemble was: -5.648% in the L15 assembly. Therefore, comparisons of the axial detector power signals were performed for the B13 and L15 assemblies. The JMCT results were located between MC21 and the measured values (see Fig. 15 Fig. 15 ). Figure 16 Figure 16 shows a comparison of the radial pin power in the axial plane between the MC21 and JMCT. The standard deviation of 99% of the fuel pin power was less than 1% (see Fig. 16(c) Fig. 16(c) ). Table 1 Table 1 shows the k eff comparison of JMCT, OpenMC, and MC21 at the different locations of the control rods and boron concentrations. Table 2 Table 2 lists the reactivity worth of the control rods at 560 0 F. Table 3 Table 3 presents a comparison of the temperature coefficients. The integrated axial power distribution of the experiment had a slight asymmetry. We guess that this was caused by the instrument tubes. Generally, the comparison shows good agreement among the JMCT, MC21, and OpenMC models [ 33 33 ] . 1 1 Critical boron concentrations (pcm) 1 1   1 1 JMCT (95% confidence) 1 1 OpenMC (95% confidence) 1 1 MC21 (95% confidence) 1 1 ARO 1 1 975 1 1 1.000479±0.00003 1 1 0.99920±0.00004 1 1 0.9992614±0.00004 1 1 D in 1 1 902 1 1 1.002174±0.00003 1 1 1.00080±0.00004 1 1 C, D in 1 1 810 1 1 1.001419±0.00003 1 1 1.00023±0.00005 1 1 A, B, C, D in 1 1 686 1 1 0.999917±0.00003 1 1 0.99884±0.00004 1 1 A, B, C, D, SE, SD, SC in 1 1 508 1 1 0.998381±0.00003 1 1 0.99725±0.00004 1 1   1 1 Control rod worth (pcm) 1 1   1 1 Measure 1 1 MC21 1 1 OpenMC 1 1 JMCT 1 1 D 1 1 938.5 1 1 788 1 1 773 1 1 771±6 1 1 770±6 1 1 C with D in 1 1 856 1 1 1203 1 1 1260 1 1 1234±7 1 1 1258±6 1 1 B with D, C in 1 1 748 1 1 1171 1 1 1172 1 1 1197±7 1 1 1162±6 1 1 A with D, C, B in 1 1 748 1 1 548 1 1 574 1 1 556±6 1 1 578±6 1 1 SE with D, C, B, A in 1 1 597 1 1 461 1 1 544 1 1 501±6 1 1 543±6 1 1 SD with D, C, B, A, SE in 1 1 597 1 1 772 1 1 786 1 1 844±6 1 1 781±6 1 1 SC with D, C, B, A, SE, SD in 1 1 597 1 1 1099 1 1 1122 1 1 1049±6 1 1 1107±6 Note: temperature coefficient under 560 0 F, T 2 = 570 0 F, T 1 = 550 0 F. 1 1 HZP criticalboron evaluation 1 1 Boron concentration 1 1 Measured(pcm/ 0 F) 1 1 MC21(pcm/ 0 F) 1 1 JMCT(pcm/ 0 F) 1 1 ARO 1 1 975 1 1 -1.75 1 1 -2.7 1 1 -2.21 1 1 D in 1 1 902 1 1 -2.75 1 1 -4.3 1 1 -2.99 1 1 C,D in 1 1 810 1 1 -8.01 1 1 -9.1 1 1 -6.22 VERA model 2 HZP condition 3 As a marked achievement of the CASL plan [ 34 34 ] , VERA core physics benchmark progression problem specifications were published on March 29, 2013 [ 5 5 ] . Figure 17 Figure 17 shows the VERA modeling by the JLAMT. Table 4 Table 4 compares the k eff and control rod worth between the JMCT and KENO-VI [ 10 10 ] . Figure 18 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. 1 1 Test result 1 1 KENO-VI 1 1 JMCT 1 1 JMCT-KENO 1 1 k -effective of 5#-2 1 1 1.000321 1 1 1.000305 1 1 -1.6 pcm 1 1 SA control rod worth (pcm) 1 1 447 ± 2 1 1 439 ± 2 1 1 -8 1 1 A control rod worth (pcm) 1 1 898 ± 2 1 1 903 ± 2 1 1 5 1 1 Boron worth (pcm/ppm) 1 1 -10.21 ± 0.02 1 1 -10.21 ± 0.02 1 1 0 1 1 Temperature coefficient (pcm/F) 1 1 -3.19 ± 0.04 1 1 -3.26 ± 0.04 1 1 -0.07 1 1 Critical boron concentration in full power (ppm) 1 1 850.5 1 1 840.5 1 1 -10 MAX=0.85%RMS=0.36% HFP condition 3 Figure 19 Figure 19 shows a comparison of the radial integral assembly power distribution between the JMCT and MC21. The maximum difference was 2.27%. AP1000 [ 6 6 ] 2 Figure 20 Figure 20 shows the CAD modeling using JLAMT. The start of the physical experiment of the first core in units A(b) and B(b) were simulated by JMCT, NECP-X [ 35 35 ] and SCAP [ 36 36 ] . Table 5 Table 5 lists the deviations of the simulation results relative to the experimental results for unit A(b). Figure 21 Figure 21 shows the deviation of the control rod worth of units A(b) and B(b). 1 1 Content 1 1 JMCT deviations 1 1 NECP-X deviations 1 1 SCAP deviations 1 1 Gray bank worth 1 1 5.7×10 -5 1 1 7.3×10 -5 1 1 -10.3×10 -5 1 1 Black bank worth 1 1 0.8×10 -5 1 1 -1.7×10 -5 1 1 3.2×10 -5 1 1 ITC/℃ 1 1 1.57×10 -5 1 1 0.28×10 -5 1 1 -0.597×10 -5 CAP1400 [ 7 7 ] 2 Reactor core 3 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 Table 6 lists the deviations of NECP-X and SCAP relative to JMCT [ 37 37 ] . Note: critical boron concentration and control bank worth of MA~MD use absolute deviations, while the others use relative deviation. 1 1 Item 1 1 JMCT results 1 1 NECP-X errors 1 1 SCAP errors 1 1 Critical boron concentration/10 -6 1 1 1224.3 1 1 -23.6×10 -6 1 1 -4.9×10 -6 1 1 Control bank worth of MA/10 -5 1 1 205.6 1 1 -5.2×10 -5 1 1 -7.6×10 -5 1 1 Control bank worth of MB/10 -5 1 1 172.8 1 1 4.0×10 -5 1 1 2.8×10 -5 1 1 Control bank worth of MV/10 -5 1 1 210.9 1 1 0.4×10 -5 1 1 -2.0×10 -5 1 1 Control bank worth of MD/10 -5 1 1 187.5 1 1 -4.8×10 -5 1 1 -6.9×10 -5 1 1 Control bank worth of M1/10 -5 1 1 498.1 1 1 -0.24% 1 1 -0.88% 1 1 Control bank worth of M2/10 -5 1 1 740.7 1 1 0.68% 1 1 0.07% 1 1 Control bank worth of AO/10 -5 1 1 1619.9 1 1 -0.14% 1 1 -1.21% 1 1 Control bank worth of SD1/10 -5 1 1 785.6 1 1 -0.87% 1 1 -1.35% 1 1 Control bank worth of SD2/10 -5 1 1 939.2 1 1 -0.05% 1 1 -0.31% 1 1 Control bank worth of SD3/10 -5 1 1 891.2 1 1 0.39% 1 1 -0.19% 1 1 Control bank worth of SD4/10 -5 1 1 788.5 1 1 -2.50% 1 1 -2.79% 1 1 Control bank worth of SD5/10 -5 1 1 463.6 1 1 2.14% 1 1 1.16% 1 1 Control bank worth of SD6/10 -5 1 1 249.8 1 1 2.68% 1 1 1.98% 1 1 Total worth of block control banks/10 -5 1 1 6976.6 1 1 -0.08% 1 1 -0.72% Shielding 3 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 Figure 22 shows the neutron dose distribution and secondary photon dose distribution of the full factory. Summary 1 A high-fidelity general-purpose 3-D Monte Carlo particle transport program, JMCT3.0, was developed for integrated simulation of nuclear systems. When the geometry zones break through ten million, the depletion regions break through millions, and parallel processors break through one hundred thousand. Advanced computer technologies, automatic CAD modeling, and visualization make the program with a high performance. The JMCT has carried out a number of simulations in reactor full-core, radiation shielding [ 38 38 ] , nuclear detection [ 39 39 ] and nuclear medicine [ 40 40 ] . A strong simulation capability was shown. In recent years, JMCT has completed the simulations of many large and complex nuclear devices, such as CFR600, ACP100, Hualong, Qingshan-I, VVER, CFETR, and HTR-10 (see Fig. 23 Fig. 23 ). The necessary technology supports the optimization of nuclear device designs. Future efforts will be forced to enhance computing efficiency. Depletion, complicating the uncertainty quantification and propagation of errors is an essential area to consider in the future. Furthermore, new algorithms need to be developed to reduce computing fees. JMCT is still evolving toward this goal, and all algorithms are being actively developed at present. References B.Y. Zhang , G. Li , L. Deng et al. , JCOGIN: A Parallel Programming Infrastructure for Monte Carlo Particle Transport . PHYSOR 2014 , Kyoto, Japan , September 28-October 3, 2014 . JCOGIN: A Parallel Programming Infrastructure for Monte Carlo Particle Transport Z. Mo , A. Zhang , X. Cao , et al ., JASMIN: a parallel software infrastructure for scientific computing . Front. Compute. Sci. China 4 ( 4 ), 480 - 488 ( 2020 ). doi: 10.1007/s11704-010-0120-5 http://doi.org/10.1007/s11704-010-0120-5 JASMIN: a parallel software infrastructure for scientific computing L. Deng , T. Ye , G. Li et al. , 3-D Monte Carlo neutron-photon transport code JMCT and its algorithms . PHYSOR 2014 , Kyoto, Japan , September 28-October 3, 2014 . 3-D Monte Carlo neutron-photon transport code JMCT and its algorithms N. Horelik , B. Herman , Benchmark for evaluation and validation of reactor simulations (BEAVRS) . RELEASE rev. 1.0.1 , MIT Computational Reactor Physics Group , July 7, 2013 . Benchmark for evaluation and validation of reactor simulations (BEAVRS) Andrew T and Godfrey , VERA Core Physics Benchmark Progression Problem Specifications. Revision 2 , CASL-U-2013-0131-002 , March 29, 2013 . VERA Core Physics Benchmark Progression Problem Specifications. F. Fausto , Zero power physics test simulations for the AP1000 PWR. CASL-U-2014-0012-001 , 2014 . Zero power physics test simulations for the AP1000 PWR. CASL-U-2014-0012-001 M.G. Zheng , J.Q. Yan , Large power advanced passive PWR CAP1400 . Shanghai: Shanghai Jiaotong University Press , 2018 . Large power advanced passive PWR CAP1400 D.J. Kelly , T.M. Sutton , S.C. Wilson , MC21 Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark Problem . Proc. Advances in Reactor Physics - Linking Research, Industry, and Education (PHYSOR 2012) , Knoxville, Tennessee , American Nuclear Society , April 15-20, 2012. MC21 Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark Problem P.K. Romano , N.E. Horelik et al ., OpenMC: A state-of-the-art Monte Carlo code for research and development . Ann. Nucl. Energy 82 , 90 - 97 ( 2015 ). doi: 10.1016/j.anucene.2014.07.048 http://doi.org/10.1016/j.anucene.2014.07.048 OpenMC: A state-of-the-art Monte Carlo code for research and development M. Stephen , KENO-VI Primer: A Primer for criticality calculations with SCALE/KENO-VI using GeeWiz . Oak Ridge National Laboratory , 2008 . KENO-VI Primer: A Primer for criticality calculations with SCALE/KENO-VI using GeeWiz R.E. MacFarlane , D.W. Muir , The NJOY Nuclear Data Processing System-Version 91. LA-12740-1994 . The NJOY Nuclear Data Processing System-Version 91. LA-12740-1994 C. Yang , T.P. Cheng , L. Deng et al ., Development of 3-D parallel first-collision source method for discrete ordinate code JSNT-S . Ann. Nucl. Energy 135 , 106942 ( 2020 ). doi: 10.1016/j.anucene.2019.106942 http://doi.org/10.1016/j.anucene.2019.106942 Development of 3-D parallel first-collision source method for discrete ordinate code JSNT-S Y.G. Fu , L. Deng , G. Li , Preliminary study on numerical solver of inhomogeneous burnup equations . Acta Phys. Sin . 67 ( 17 ), 172802 ( 2018 ). doi: 10.7498/aps.67.20172650 http://doi.org/10.7498/aps.67.20172650 ( in Chinese ) Preliminary study on numerical solver of inhomogeneous burnup equations Y. Ma , Y.G. Fu , G.M. Qin , The design of JLAMT: An aided tool for large-scale complex physical modeling . Advance in Intelligent Systems and Computing , 809 , 877 - 883 ( 2019 ). doi: 10.1007/978-3-319-95588-9_74 http://doi.org/10.1007/978-3-319-95588-9_74 The design of JLAMT: An aided tool for large-scale complex physical modeling Y. Cao , Z.Y. Mo , L. Xiao et al ., Efficient visualization of high-resolution virtual nuclear reactor . J. Visualization 21 , 857 - 871 ( 2018 ). doi: 10.1007/s12650-018-0487-1 http://doi.org/10.1007/s12650-018-0487-1 Efficient visualization of high-resolution virtual nuclear reactor X-5 Monte Carlo team , MCNP-A General Monte Carlo Code for N-Particle Transport Code . version 5 manual, LA-UR-03-1987 , April 24, 2003 . MCNP-A General Monte Carlo Code for N-Particle Transport Code HEMANM , TRKOV A , ENDF-6 formats manual, data formats and procedures for the evaluated nuclear data file ENDF/B-VI and ENDF-VII. BNL , 2009 . ENDF-6 formats manual, data formats and procedures for the evaluated nuclear data file ENDF/B-VI and ENDF-VII. BNL Z.G. Ge , R.R. Xu , H.C. Wu et al. , CENDL-3.2: The new version of Chinese general purpose evaluated nuclear data library . EPJ Web of Conferences 239 , 09001 ( 2020 ). doi: 10.1051/epjconf/202023909001 http://doi.org/10.1051/epjconf/202023909001 CENDL-3.2: The new version of Chinese general purpose evaluated nuclear data library R.W. Roussin , BUGLE-80 Coupled 47 Neutron and 20 Gamma-Ray, P 3 Cross Section Library for LWR Shielding Calculations. DLC-75 ( 1980 ). BUGLE-80 Coupled 47 Neutron and 20 Gamma-Ray, P3 Cross Section Library for LWR Shielding Calculations. DLC-75 L.L. Carter , E.D. Cashwell , Particle Transport Simulation with the Monte Carlo Method, ERDA Critical Review Series. TID-26607 ( 1975 ). Particle Transport Simulation with the Monte Carlo Method, ERDA Critical Review Series. TID-26607 G. Li , B. Zhang , L. Deng et al. , Combinatorial Geometry Domain Decomposition Strategies for Monte Carlo Simulations . M&C 2013 , Sun Valley, Idaho , May 5-9 , pp. 434 - 443 ( 2013 ). Combinatorial Geometry Domain Decomposition Strategies for Monte Carlo Simulations R. Procassini , D. Cullen , G. Greenman et al. , New capabilities in MERCURY: a modern Monte Carlo particle transport code . M&C+SNA , Monterey California , April 15-19, 2007 . New capabilities in MERCURY: a modern Monte Carlo particle transport code D.J. Kelly , T.M. Sutton , S.C. Wilson , MC21 Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark Problem . Proc. Advances in Reactor Physics—Linking Research, Industry, and Education (PHYSOR 2012) , Knoxville, Tennessee , April 15-20, 2012 . MC21 Analysis of the Nuclear Energy Agency Monte Carlo Performance Benchmark Problem D.H. Shangguan , G. Li , B.Y. Zhang et al ., Uniform tally density-based strategy for efficient global tallying in Monte Carlo criticality calculation . Nucl. Sci. Eng . 182 ( 4 ), 555 - 562 ( 2016 ). doi: 10.13182/NSE15-32 http://doi.org/10.13182/NSE15-32 Uniform tally density-based strategy for efficient global tallying in Monte Carlo criticality calculation X.G. Liu , L. Deng , Z.H. Hu et al ., Study of on-the-fly Doppler broadening in JMCT Program . Acta Phys. Sin . 65 ( 9 ), 092501 ( 2016 ). doi: 10.7498/aps.65.092501 http://doi.org/10.7498/aps.65.092501 ( in Chinese ) Study of on-the-fly Doppler broadening in JMCT Program R. Li , L.Y. Zhang , D. Shi et al ., Criticality search of soluble boron iteration in Monte Carlo code JMCT . Energy Procedia , 127 , 329 - 334 ( 2017 ). doi: 10.1016/j.egypro.2017.08.118 http://doi.org/10.1016/j.egypro.2017.08.118 Criticality search of soluble boron iteration in Monte Carlo code JMCT Z. Zheng , Q.L. Mei , L. Deng , Study on variance reduction technique based on adjoint Discrete Ordinate method . Ann. Nucl. Energy 112 , 374 - 386 ( 2018 ). doi: 10.1016/j.anucene.2017.10.028 http://doi.org/10.1016/j.anucene.2017.10.028 Study on variance reduction technique based on adjoint Discrete Ordinate method Z. Zheng , M.Q. Wang , L. Deng , Application of a 3D Discrete Ordinates-Monte Carlo coupling method to deep-penetration shielding calculation . Nucl. Eng. Des . 326 , 87 - 96 ( 2018 ). doi: 10.1016/j.nucengdes.2017.11.005 http://doi.org/10.1016/j.nucengdes.2017.11.005 Application of a 3D Discrete Ordinates-Monte Carlo coupling method to deep-penetration shielding calculation L. Deng , Z.H. Hu , R. Li et al ., The coupled neutron transport calculation of Monte Carlo multi-group and continuous cross section . Ann. Nucl. Energy 127 , 433 - 436 ( 2019 ). doi: 10.1016/j.anucene.2018.12.032 http://doi.org/10.1016/j.anucene.2018.12.032 The coupled neutron transport calculation of Monte Carlo multi-group and continuous cross section M.B. Emmett , MORSE-CGA - A Monte Carlo radiation transport code with array geometry capability. ORNL-6174 , 1985 . MORSE-CGA - A Monte Carlo radiation transport code with array geometry capability. ORNL-6174 D. She , B. Xia , J. Guo et al ., Prediction calculations for the first criticality of the HTR-PM using the PANGU code . Nucl. Sci. Tech . 32 , 90 ( 2021 ). doi: 10.1007/s41365-021-00936-5 http://doi.org/10.1007/s41365-021-00936-5 Prediction calculations for the first criticality of the HTR-PM using the PANGU code D.J. Kelly , B.R. Herman et al. , Analysis of select BEAVRS PWR benchmark cycle 1 results using MC21 and OpenMC . PHYSOR 2014 , Kyoto, Japan , September 28-October 3, 2014 . Analysis of select BEAVRS PWR benchmark cycle 1 results using MC21 and OpenMC G. Li , L. Deng , B.Y. Zhang et al ., JMCT Monte Carlo analysis of BEAVRS benchmark: hot zero power results . Acta Phys. Sin . 65 ( 5 ), 052801 ( 2016 ). doi: 10.7498/aps.65.052801 http://doi.org/10.7498/aps.65.052801 ( in Chinese ) JMCT Monte Carlo analysis of BEAVRS benchmark: hot zero power results P.J. Turinsky , D.B. Kothe , Modeling and simulation challenges pursued by the Consortium for Advanced Simulation of Light Water Reactors (CASL) . J. Computational Phys . 313 , 367 - 376 ( 2016 ). doi: 10.1016/j.jcp.2016.02.043 http://doi.org/10.1016/j.jcp.2016.02.043 Modeling and simulation challenges pursued by the Consortium for Advanced Simulation of Light Water Reactors (CASL) J. Chen , Z.Y. Liu , C. Zhao et al ., A new high-fidelity neutronics code NECP-X . Ann. Nucl. Energy 116 , 417 - 428 ( 2018 ). doi: 10.1016/j.anucene.2018.02.049 http://doi.org/10.1016/j.anucene.2018.02.049 A new high-fidelity neutronics code NECP-X L.H. Peng , C.T. Tang , W.Y. Yang , Development of advanced neutronics code SCAP-N for reactor core high fidelity simulation . Nuclear Power Engineering 42 ( 2 ), 213 - 218 ( 2021 ). ( in Chinese ) Development of advanced neutronics code SCAP-N for reactor core high fidelity simulation L.H. Peng , B. Yang , C.T. Tang et al ., Zero power physics test high fidelity simulation for first core of Guo He One (CAP1400) reactor . High Power Laser and Particle Beams 34 ( 2 ), 026002 ( 2022 ). doi: 10.11884/HPLPB202234.210372 http://doi.org/10.11884/HPLPB202234.210372 ( in Chinese ) Zero power physics test high fidelity simulation for first core of Guo He One (CAP1400) reactor L. Deng , R. Li , Q.X. Din et al ., Qinshan-I reactor shielding simulation and sensitivity analysis based on JMCT monte carlo code . Nuclear Power Engineering 42 ( 2 ), 173 - 179 ( 2021 ). ( in Chinese ) Qinshan-I reactor shielding simulation and sensitivity analysis based on JMCT monte carlo code L. Deng , R. Li , X. Wang et al ., Monte Carlo simulation technology based on characteristic γ-ray spectrum analysis . Acta Phys. Sin . 69 ( 11 ), 112801 ( 2020 ). doi: 10.7498/aps.69.20200279 http://doi.org/10.7498/aps.69.20200279 ( in Chinese ) Monte Carlo simulation technology based on characteristic γ-ray spectrum analysis G. Li , L. Deng , Optimized voxel model construction and simulation research in BNCT . Chinese Phy. C 30 ( 2 ), 171 - 177 ( 2006 ). Optimized voxel model construction and simulation research in BNCT All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Li Deng, Ze-Yao Mo, Gang Li, Bao-Yin Zhang, Rui Li, Yuan-Gang Fu, Dun-Fu Shi, Xin Wang, Dan-Hu Shangguan, Ling-Yu Zhang, Peng Liu, Ze-Hua Hu, Yan Ma, Sheng-Cheng Zhou and Jing-Wen Shen. The first draft of the manuscript was written by Li Deng and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript. 10.1007/s41365-022-01092-0 This work was supported by the National Natural Science Foundation of China (Nos.11805017 and 12001050). 2022-04-28 2022-07-30 2022-08-01 2022-08-27 2022-09-26 © China Science Publishing & Media Ltd. (Science Press), Shanghai Institute of Applied Physics, the Chinese Academy of Sciences, ChineseNuclear Society and Springer Nature Singapore Pte Ltd. 2022 This is a post-peer-review, pre-copyedit version of an article published in Nuclear Science and Techniques. The final authenticated version is available online at: http://dx.doi.org/10.1007/s41365-022-01092-0 http://dx.doi.org/10.1007/s41365-022-01092-0 . 2022 Nuclear Science and Techniques 08 33