Progress and prospects of calculation methods for radiation shielding

Abstract

Progress in calculation methods for radiation shielding are reviewed based on the activities of research committees related to radiation shielding fields established in the Atomic Energy Society of Japan. A technological roadmap for the field of radiation shielding; progress and prospects for specific shielding calculation methods such as the Monte Carlo, discrete ordinate Sn transport, and simplified methods; and shielding experiments used to validate calculation methods are presented in this paper.

Keywords:

• radiation shielding
• gamma-ray
• neutron
• Monte Carlo
• Sn transport method
• simplified calculation method
• benchmark

1. Introduction

Radiation shielding is indispensable for the safety of nuclear reactors and spent-fuel casks, fuel-cycle, radioisotope-handling, accelerator and nuclear-fusion-test facilities, and so on, where radiation and radioactive materials are used. These facilities and systems are used for industry, research, medical, treatment, and diagnostics. Safety regulations for these systems and facilities are essential to ensure the health and safety of the public. In 1964, the radiation shielding research committee was established within the Atomic Energy Society of Japan (AESJ) to collect data and analyze information regarding radiation sources and behaviors in various radiation facilities. Since its establishment, the safety and reliability of radiation shielding in these facilities has dramatically improved. In accordance with the recent advances in research and the development of computer technologies in each field, shielding analyses and design methods have also become sophisticated, and the applicability of high energy accelerators to new facilities has been discussed. And, the use of radioactivity is becoming more popular in medical and industrial fields, the range of areas to which radioactivity calculations and relevant data apply is increasing.

In this review, the “Technological Roadmap in the Radiation Shielding Field”, prepared by the “Radiation Shielding research committee established in the AESJ” is first considering future prospects, technological levels need to be confirmed in light of the current needs and remaining issues for the radiation shielding technological development need to be sorted out. Subsequently, the progress and prospects of individual calculation methods for radiation shielding, such as the Monte Carlo, discrete ordinate Sn transport, and simplified calculation methods, are presented. Finally, shielding experiments used to establish benchmarks for shielding calculation methods are presented.

2. Technological roadmap in radiation shielding field

2.1. Outlines of the roadmap in radiation shielding

In consideration of the current applied fields and working status, the “Radiation Shielding” research committee classified the current issues and future trends into the following categories: reactor facilities, reprocessing/fabrication facilities, decommissioning/low-level waste, transport/storage facilities, nuclear fusion reactors, accelerators, medical equipment, aerospace, industry utilization, and others. The committee compiled the technological issues to be tackled in the future in these applied fields and facilities.

Four years were needed to chart out this technological roadmap. However, substantial changes were made in the final phase of the roadmap owing to the occurrence of the Great East Japan Earthquake and the related accidents of Tokyo Electric Power Company's Fukushima No. 1 nuclear power plant (hereinafter referred to as the Fukushima No. 1 plant). The contents in Sections 2.2–2.4 were considered before the accident of the Fukushima No. 1 plant. Section 2.5 summarizes the situational changes for individual fields after the Fukushima No. 1 plant accident.

2.2. Current issues and future prospects regarding radiation shielding

Current issues and future prospects regarding radiation shielding are classified into the following categories: reactor facilities, reprocessing/fabrication facilities, decommissioning/low-level waste, transport/storage facilities, nuclear fusion reactors, accelerators, medical equipment, aerospace, industry utilization, and others. The common issues and prospects among these are summarized below.

In association with the progress and improvement of radiation shielding analysis technologies, the uncertainties attributed to the analysis have decreased. Use of the latest insights and knowledge can contribute to more rational and economic shielding designs. However, concerns about “radiation” remain among the general public; furthermore, determining the appropriate margin settings for the shielding design is difficult because insights into the engineering of radiation behavior are imperfect. Currently, the measures taken to introduce the latest insights into the safety regulations are inadequate. Overcoming this situation requires analysis standardization and validation of analysis methods, accumulation of data, and maintenance of irradiation facilities. Moreover, international cooperation and information exchange systems must be enhanced. To achieve these objectives, the development of globally compatible Japanese analysis codes, nuclear data, and enhancement of technology improvement are indispensable.

A common challenge across all fields is ensuring that the accumulated radiation shielding technologies and insights are passed on to future generations. Standardization and development of guidelines for the shielding evaluation method are necessary for training and developing the abilities of young researchers and engineers, and for transmitting the knowledge of shielding calculation codes. Additional tasks to prevent the dissipation of accumulated insights and knowledge about shielding evaluation technologies and goal setting to promote sequential basic research program are necessary.

Furthermore, in consideration of the special characteristics of licensed reactors, high reliability and rationality should be assured in the design method and process, and their evaluation method.

The research committee examined issues by focusing on the field of radiation shielding. Other issues related to the fields intimately connected with radiation shielding such as radiation monitoring, irradiation level evaluation, and radiation exposure doses remain to be considered. For example, when the International Commission on Radiological Protection (ICRP) published new recommendations in 2007 [1], a part of both the radiation weighting factor and tissue weighting factor, which are necessary for the calculation of the effective dose, were changed. In 2012, coefficients for conversion to the effective dose reflecting this change were provided [2].

Once the new recommendations are adopted in the national law and ordinance, effective dose rate is calculated by multiplying the dose conversion coefficient to the spectral data in the detailed calculation method. In the simplified calculation method, which uses shielding calculation constants that depend on the dose conversion coefficient, new shielding calculation constants like gamma ray buildup factors must be prepared.

As stated above, concrete action plans for creating the technological roadmap must be provided to solve the aforementioned issues while identifying other issues by broadening the scope to all fields of radiation engineering.

2.3. Fact-finding study of radiation shielding field

2.3.1. Identification of issues in individual fields

To understand the types of issues in individual fields, fact-finding studies were conducted to extract issues to be submitted to the committee. The resulting issues were summarized in three different categories: “Main issues”, “Design/analysis status”, and “other status”.

This study has clarified technological issues regarding individual radiation shielding fields. These issues are used as the base data for category classification in the process of preparing a roadmap. The results of issue identification are omitted here. Although the issues in each field differ slightly in type, those regarding analysis codes and handbooks, the training of engineers, and technology transfer to successor have been identified as common issues.

2.3.2. Questionnaire results of fact-finding study

The issues facing the respective fields were understood to some extent in Section 2.3.1. To obtain more precise information, a questionnaire survey was conducted. This survey was conducted from May through June of 2010.

In the survey, 21(in person) effective answers were obtained. Figure 1 shows the summarized results in the form of pie charts. Figure 1.

Pie charts of summarized results for questionnaire survey.

The obtained tendencies are summarized as follows:

• (1) Overall

• Most radiation shielding departments consist of a small number of engineers (two to five).

• Analyses are mainly conducted inside the organization and analysts have the ability to evaluate the results.

• Personnel are mainly composed of engineers with experience.

• Although securing the necessary personnel has been possible thus far, it will be difficult in the future.

• The training of engineers and technology transfer involve the methods of actual operations and training by personnel such as OJT (On Job Training) and lectures.

• Inadequate training is due to insufficient personnel and busy work schedule.

• References are mainly based on documents prepared in Japan.

• MCNP [3,4] is used predominantly as the calculation code.

• The nuclear data JENDL [5] is widely used, and the most suitable version is selected in accordance with the issues.

• FORTRAN is used as the program language for all institutes.

• (2) Characteristics of each institute

• Personnel composition and future prospects are different for each institute.

• The characteristics of the codes to be used vary depending on the purpose – whether the codes are used for a reactor facility or for an accelerator.

2.4. Technological roadmap for radiation shielding field

2.4.1. Required items for radiation shielding field

Radiation shielding roadmaps were reviewed on the basis of identifying the issues of each field and the fact-finding study. The original plan was to describe the required items of the radiation shielding field in connection with the roadmaps of other fields. However, owing to the occurrence of the Great East Japan Earthquake on 11 March 2011, future prospects for the nuclear field as a whole have become vague. Consequently, mutually linking the roadmaps of various fields has become difficult.

Therefore, selected items from past reviews are summarized by classifying them into the following three categories:

1. Items related to standardization (verification and validation [V&V])

2. Items related to improvement of evaluation accuracy

3. Items related to technological hand-down

To prepare the roadmap for the radiation shielding field, each item is addressed as “short-term (5 years)”, “mid-term (5 to 10 years)”, or “long-term (10 to 15 years)”. Selected items from the questionnaire are thus described with respective terms in the roadmap.

2.4.2. Status quo of each item

(1) Items related to standardization (V&V)

• Standard data for the gamma ray buildup factor are now being arranged by AESJ. In these data, photon cross section data and the buildup factor calculation code are unified; the shielding thickness is expanded from 40 mean-free paths (mfp) to 100 mfp; and 5 different doses (exposure dose, absorbed dose, effective dose for antero-posterior (AP) irradiation, effective dose for isotropic (ISO) irradiation, and 1 cm dose equivalent) are considered as the buildup coefficient. Standardization of the typical elemental composition of the shielding materials will be considered. In this regard, a standardization committee of the American Nuclear Society (ANS) has started reviewing standard data for the photon attenuation coefficient and the gamma ray buildup coefficient, both of which were established for engineering materials in 1991 [6].

• With regard to the Monte Carlo calculation code used for the shielding calculations in licensing evaluations, guidelines to use the code for shielding calculation of cask are considered [7].

(2) Items related to improvement of evaluation accuracy

• The particle and heavy iron transport calculation code PHITS2 [8] and the cross sectional data set of neutrons and photons based on JENDL-4 library have been jointly developed by the Japan Atomic Energy Agency (JAEA), High Energy Accelerator Research Organization (KEK), and Research Organization for Information Science and Technology (RIST), and they are now available. Consequently, the number of users in the medical application field has increased.

• In the evaluation of the effective-dose conversion coefficient for radioisotopes, which are given in the Radioisotope Notebook issued by the Japan Radioisotope Association the lower energy limit of the considering gamma-ray and X-ray emission was lowered from 30 keV to 10 keV [9]. In accordance with this lowering, the “Shielding Calculation Operational Manual for Radiation facilities”, issued by the Nuclear Safety Technical Center, is being revised [10].

• Although the values described in the TecDoc Report of the IAEA [11] are most commonly used in the dose evaluation at accidents, however, in case of some data, the details and processes by which they are derived are unspecified. Therefore, the data, for example, used for decontamination evaluation, must be re-obtained with the latest calculation method in dose evaluation. For example, when the representative nuclides are uniformly distributed on the ground, the scattered radiation obtained by use of the point kernel integral calculation code, QAD [12], is approximately double the comparable obtained by the Monte Carlo calculation code EGS [13] or PHITS [8].

(3) Items related to technology transfer to the successors

• Although the neutron shielding handbook prepared by the radiation shielding related research committees and the gamma-ray shielding handbook [14,15] were issued 20 years ago, they are still used even today. Thus, these handbooks need to be revised to introduce the latest information.

2.4.3. Roadmap for radiation shielding field

Table 1 presents. a roadmap of the radiation shielding field. For (1)items related to standardization (V&V) and (2) improvement of evaluation accuracy, the described contents are mainly about standardization or higher accuracy of the analysis code/library. For (3) items related to technological hand-down, the contents are concerned with the handbooks, training courses, and maintenance of the test facilities.

Table 1. Technological roadmap of radiation shielding field.

In particular, as shown in Section 2.4.2, Japanese code development and shielding handbook revision need to be tackled with high priority.

2.5. Status of respective fields after Fukushima No. 1 nuclear power plant accident

2.5.1. Reactor, fuel field

The accident of Fukushima No.1 plant has significantly affected the field of nuclear energy as a whole, including reviews not only of light water reactors but also of fuel cycles. The status of nuclear power plants in Japan is that they will all be shut down by April 2012 for periodic inspection. With regard to the resumption of these plants, numerous steps and time are required for preparing emergency safety measures against tsunami and earthquakes, examining and reporting stress tests, and gaining permission from local governments.

Currently, reviews of the Basic Energy Plan, which includes nuclear power generation, nuclear safety regulations and systems, and nuclear disaster prevention systems, are promoted; in addition, the safety measures against severe accidents at light water reactors are being strengthened. In this regard, radical and drastic changes are expected.

Systematic removal of radioactive cesium across a broad region is promoted to allow residents to return home in Fukushima. In addition, detailed dose monitoring, reviews of restricted zones, and various support policies after the accident are promoted for this purpose. To achieve decommissioning of the Fukushima No. 1 plant, safety assurance measures for the public and radiation workers (reactor stabilization, on-site decommissioning and decontamination, study of relevant technologies) are now underway in accordance with the “Basic concept of mid-term safety assurance [16]”.

In light of the above situation, the licensing and building of new nuclear power plants and fuel cycle facilities are difficult in Japan. In the meantime, except for international projects, the required work will mainly consist decommissioning of nuclear plants that are slated to be abolished, building intermediate storage facilities for the removed soil of decontamination, and disposing of or storing debris and fuels related to the Fukushima No. 1 plant.

Because these various impacts of the Fukushima No. 1 plant accident and nuclear safety issues draw significant public attention, the importance of safety and its evaluation technology in relation to severe nuclear reactor accidents is rising. Radiation shielding technologies are also of great importance for environmental dose evaluation during severe accidents.

Stated simply, excepting international projects, the expectations for new technologies, new data, and increased demand for training new personnel are low for shielding in the reactor and fuel field. For the time being, the status quo must be maintained. However, operations regarding radiation and radioactive materials, including the above-described nuclear safety evaluations, are rather increasing owing to their wide range of applicability. From this viewpoint, the role of the radiation shielding and radiation exposure field has increased.

2.5.2. Waste/transport/storage field

Owing to the Fukushima No. 1 plant accident, the following items have emerged as the new issues in the field of waste/transport/storage.

• Storage and transport of radioactive waste generated by decontamination of vast areas of soil contaminated with radioactive materials that were released to the fallout during the accident.

• Storage and transport of radioactive waste such as debris containing high level radioactivity in restricted zones.

• Storage and transport of radioactive waste generated by decommissioning of contaminated reactor facilities.

• Transport of damaged fuels in reactors and in spent fuel storage pools.

Decontamination operations have already commenced in the restricted zones. The radioactive waste is first collected in a certain place, and transported to intermediate storage facilities. Mid-term or long-term operations are expected for dealing with debris in the restricted zones. All radioactive wastes generated in this accident differ from the conventional assumptions, and thus evaluation of their sources is a challenging work. However, the situations can be sufficiently managed and addressed using the existing shielding technologies.

Whereas decontamination operations are generally implemented for a certain range of areas, the dose contributions from non-decontaminated zones cannot be ignored. To address this problem, multiple institutes such as JAEA and universities provide systems that can evaluate dose reduction effects and easily calculate the dose after decontamination by using a matrix in which the dose contribution of the zones away from the evaluation point is calculated in advance.

Transport of the damaged fuel is a mid- and long-term issue because it is to be implemented after the technical issues regarding taking out of the molten core debris are solved. As the example of TMI-2 showed, certain issues of identifying the radioactivity inventory must be addressed in the area. Apparently, designs for containment, etc., will be first needed. In developing these designs, methods for identifying of the radioactivity inventory will likely need to be discussed in terms of where to set the capacity limit of containment. Conventional shielding evaluation technologies are able to containment safety evaluation.

2.5.3. Nuclear fusion field

In the field of nuclear fusion, the impacts of the Fukushima No. 1 plant accident have been minor. With their broad approaches, the International Thermonuclear Experimental Reactor (ITER) project [17] and the BA project [18] continue to make steady progress. However, nuclear fusion reactors and their safety also tend to draw significant attention. These safety issues are becoming the prominent topics of discussion at academic organizations and meetings on nuclear fusion. For nuclear fusion reactors, an issue of increasing significance, in addition to general radioactive shielding, is the confinement of tritium used as fuels. Radiation shielding researchers and tritium researchers will need to work together more closely in the future.

2.5.4. Accelerator, radiation application, and other fields

The impacts of the Fukushima No. 1 plant accident have been minor in other fields such as those of accelerators and radiation applications. In relation to the earthquake disaster, a plan to build synchrotron radiation facilities in northeast Japan is under consideration.

2.5.5. Summary of each field

The Fukushima No.1 plant accident crucially affected the fields related to nuclear technology in Japan, with varying degrees of impact depending on the field. Currently, reviews of conventional nuclear policies are underway.

As representatives of the country responsible for the accident, Japanese experts in radiation exposure have a responsibility to address and share the lessons learned from the Fukushima No. 1 plant accident with the global community. Applying these lessons to the radiation shielding field is a significant challenge. The world's perception of nuclear power generation, which is the mainstay of nuclear energy utilization, has dramatically changed after the Great East Japan Earthquake. Under such circumstances, the hand-down of accumulated evaluation technologies regarding radiation shielding has become an important issue.

In this view, coordination of existing technologies is indispensable for sharing the lessons learned from the Fukushima No. 1 plant accident and to conducting engineer training as well as technological hand-down. To achieve these objectives, the first step is to create Japanese original shielding calculation codes and coordinate shielding handbooks.

3. Monte Carlo method

The Monte Carlo method solves the Boltzmann transport equation as a Fredholm integral equation for particle collision probability with nuclear data in low-energy regions (<20 MeV). In this method, the behavior of particles in materials is traced by using a pseudorandom number generator. Thus, a more precise name for the method would be the Markov chain Monte Carlo method. A Markov chain is a mathematical system that undergoes transitions between a finite number of possible states.

The method is capable of treating highly complex three-dimensional (3D) configurations. The continuous treatment of energy, space, and angle obviates discretization errors such as those resulting from the use of multi-group approximation for energy and of segmentation for angle. However, the method generally requires extensive computing time as compared with other methods. Because particles are traced and counted at the estimation points or regions, the fractional standard deviation of the results should be reduced to obtain precise estimation results. The performance strongly depends on the calculation speed of computers. To overcome the problem, many variance reduction techniques such as the regional importance, weight window, and forced collision methods have been developed, and efforts to develop these are ongoing.

Recently, the method has become remarkably popular because of the dramatic progress of computers. The method can be used for small-scale simulations (e.g., response calculation of detectors) in personal computers, and for large-scale simulations (e.g., of designs for nuclear facilities) in supercomputers with many parallel processors.

In concert with the progress of computers, Monte Carlo codes have been actively developed over the past 40 years. Currently, the MVP [19–21] and MCNP [3,4] codes, developed in Japan and USA for neutron and photon simulation in the energy region of less than 20 MeV, are popular worldwide. Unfortunately, MCNP is under export control by the US government, and its use and distribution are strictly limited.

The progress of computers has extended the energy region and types of particles treated in Monte Carlo codes. In the high-energy region, the Monte Carlo code is usually applied to trace particles generated in INC interactions, because the generation process is complicated and involves multiple particles. In this energy region, charged particles should be traced as well because they have long and non-negligible ranges in materials. Precise trace of their behaviors is difficult because charged particles lose their energy by ionization and thousand-fold scatterings occur. To overcome this difficulty, the continuous slowing down approximation has been applied to the codes.

The extension of energy and particle regions has made the Monte Carlo codes more popular. Consequently, their fields of application have explosively expanded from shielding design and dose estimation around nuclear facilities to accelerator designs for high-energy physics and medical facilities, dose estimation for diagnostics and treatment in the medical field, analyses of nuclear reaction mechanisms in astrophysics, and dose estimation in space craft and aircraft as well as radiation damage estimation of materials such as semi-conductors in the space environment. To fulfill the demands of the application fields, EGS5 [13] and PENELOPE [22] have been actively developed for electromagnetic cascade simulation. Additionally, PHITS, MARS [23,24], MCNP6 [4], FULKA [25,26], and GEANT [27] have been developed as comprehensive code systems.

In this section, MVP, EGS5, and PHITS are introduced as representative examples of recent Japanese work on Monte Carlo code development.

3.1. MVP/GMVP [16–18]

MVP/GMVP is a general purpose continuous-energy Monte Carlo code for neutron and photon transport developed at JAEA. MVP is a continuous-energy version, and GMVP is a multi-group version of MVP. MVP has been developed for fast and precise simulations on nuclear reactor analysis, criticality safety, and shielding design since the late 1980s. After the first version was released in 1994, the second version was extensively modified and released for worldwide use in 2005. MVP has been used especially in the field of nuclear reactor analysis and criticality safety.

The code solves eigenvalue and fixed-source problems; neutron, photon, and neutron-photon coupled transport problems; and time-dependent problems. Combinatorial geometry is employed for geometrical modeling with multiple-lattice capability. For neutron reactions, all reactions given in the ENDF-6 format are explicitly treated in the energy range of 10^–5 eV–20 MeV. Unresolved resonances are treated by the probability table method. Thermal scattering is treated with scattering law or free gas models. For photon reactions, coherent and incoherent scattering, the photoelectric effect, and pair production are explicitly treated in the energy range of 1 keV–100 MeV. Fluorescent X-rays and Bremsstrahlung photons can be also taken into account.

In the following sections, the latest unique capabilities of MVP are described.

3.1.1. Statistical geometry model

The statistical geometry model (STGM) has been developed for high-temperature gas-cooled reactors (HTR). In HTR, coated fuel particles used as the fundamental fuel elements are randomly packed in compacts or pebbles. In the case of pebble bed reactors, fuel pebbles are also randomly stacked in the core. The position of each sphere is sampled from a probability density function called the nearest neighbor distribution (NND). The NND is obtained by a hard sphere packing code MCRDF [28], either exactly or by a theoretical approximation based on a statistically uniform distribution as follows: (1) $$\begin{array}{c}\hfill \frac{\mathit{dNND}\left(\mathit{r}\right)}{dr}=\frac{3}{2}·\frac{f_p}{1-f_p}·exp\left(-\frac{3}{2}·\frac{f_p}{1-f_p}·r\right),\end{array}$$(1) where r is the distance in unit of the diameter of the fuel particles and f_p is the packing fraction.

Figure 2 shows the schematic diagram of the improved STGM. For a neutron entering a stochastic mixture region, the distance to a sphere is sampled from an NND and the distance to a collision point is also sampled. If the distance to the collision point is larger than that to the sphere, the neutron is assumed to enter the sphere. Subsequently, the type of sphere into which the neutron enters is chosen according to the existence ratio. In the case of two spheres in the region, as shown in this figure, sphere 1 is selected first, and the neutron takes a usual random walk in sphere 1. If the neutron exits sphere 1, it repeats the same sequence as before. In this case, the neutron is assumed to enter sphere 2, after which it takes a random walk and exits sphere 2. If the distance to a collision point is smaller than that to a sphere, the neutron collides with some nuclide in the region. If the distance to a sphere is larger than that to the boundary, the neutron exits from the mixture region.

Figure 2. Schematic diagram of the improved statistical geometry model.

3.1.2. Calculation at arbitrary temperature

MVP libraries are generated from evaluated nuclear data with the LICEM code system [29]. By using the LICEM system, the cross section library of MVP can be produced from these evaluated nuclear data at any temperature. The LICEM system also generates interpolation parameters between different temperatures for thermal scattering data and data in the unresolved resonance region. To avoid the computational burden of processing probability tables for the unresolved resonance region and large storage requirements for point-wise cross sections, the ART code [30] has been also developed for MVP. ART generates cross section data at an arbitrary temperature by Doppler broadening for the resolved resonance region. The interpolation between temperatures is performed for unresolved resonance and thermal scattering data. Because ART is also integrated into MVP, the only requirement is to specify a temperature for each material in the input data. MVP then automatically generates cross sections at specified temperatures before the random walk begins.

3.1.3. Burn-up calculation

To perform burn-up calculations, MVP is coupled with a depletion module based on the predictor corrector method. The system is released as the MVP-BURN code [31].

Figure 3 shows, as an example of an MVP-BURN calculation, the k-infinity value as a function of burn-up in a MOX (mixed oxide)-fueled BWR (boiling water reactor) assembly. In the calculation, a 10 × 10 lattice consists of MOX fuel pins and uranium-oxide gadolinia pins. The results of calculations with MVP-BURN and with the VMONT multi-group Monte Carlo code [32] by Hitachi are compared. Both results show good agreement less than 0.3% through burn-up.

Figure 3. k-infinity as a function of burn-up in a MOX-fueled BWR assembly of OECD/NEA/NSC benchmark on MOX-BWR.

3.2. EGS5 [13]

The EGS (Electron-Gamma Shower) code system is a general-purpose package for Monte Carlo simulations on the coupled transport of electrons and photons in the energy region from a few keV up to several hundred GeV. The first version of EGS was developed for high-energy physics at SLAC (Stanford Linear Accelerator Center). Since the release of EGS3 [33] with a manual in 1978, it became more popular to various research fields in the general sciences because EGS3 included low-energy electron-photon transport. In 1985, EGS4 [34] was released by SLAC, KEK, and NRCC (National Research Council Canada) as a comprehensive electron and photon transport code, and was widely used in various research fields such as medical physics for about 20 years. After release of EGS4, improvements continued to be made mainly for the low-energy transport region by SLAC, KEK, and Michigan University. Finally, EGS5 [13] was released in 2001.

In the EGS5 Code System, the radiation transport of electrons/positrons or photons in any element, compound, or mixture can be simulated. The data preparation package, PEGS5, creates the data to be used by EGS5, using cross section tables for elements 1 through 100 from nuclear data libraries such as PHOTOX [35]. The dynamic range of photon energies is between 1 keV and several thousand GeV. In EGS5, the geometry is specified by a user-written subroutine with auxiliary geometry routines for planes, cylinders, cones, spheres, and so on, and the MORSE-CG Combinatorial Geometry package [36]. EGS5 takes into account various physical processes: (1(1) $$\begin{array}{c}\hfill \frac{\mathit{dNND}\left(\mathit{r}\right)}{dr}=\frac{3}{2}·\frac{f_p}{1-f_p}·exp\left(-\frac{3}{2}·\frac{f_p}{1-f_p}·r\right),\end{array}$$(1) ) Bremsstrahlung production; (2(2) $$\begin{array}{ccc}& & \nabla ·\mathit{\Omega }\Phi (\mathit{r},\mathit{\Omega },E)+\sum _t(\mathit{r},E)\Phi (\mathit{r},\mathit{\Omega },E)\hfill \\ & & ={\int }_{E^{\text{'}}}d{E}^{\text{'}}{\int }_{\mathit{}{\mathbf{\Omega }}^{\text{'}}}d\mathit{}{\mathbf{\Omega }}^{\text{'}}\sum _S\mathit{r},E^{\text{'}}\to E,\mathit{}{\mathbf{\Omega }}^{\text{'}}\to \mathit{\Omega })\Phi (\mathit{r},\mathit{}{\mathbf{\Omega }}^{\text{'}},E^{\text{'}})\hfill \\ & & +S(\mathit{r},\mathit{\Omega },E)\hfill \end{array}$$(2) ) positron annihilation in flight and at rest; (3(3) $$\begin{array}{ccc}& & {\Phi }_L\mu dy+{\Phi }_B\eta dx-{\Phi }_R\mu dy\hfill \\ & & -{\Phi }_T\eta dx-{\Sigma }_t{\Phi }_{C}dv=Q_{s}dv+Sdv\hfill \end{array}$$(3) ) Moliére multiple scattering (Coulomb scattering from nuclei); (4(4) $$\begin{array}{c}{\Phi }_C=\alpha {\Phi }_R+(1-a){\Phi }_{L}0.5\leqq \alpha \leqq 1\hfill \\ {\Phi }_C=\beta {\Phi }_B+(1-\beta ){\Phi }_{T}0.5\leqq \beta \leqq 1\hfill \end{array}$$(4) ) Møller (e⁻e⁻) and Bhabha (e⁺e⁻) scattering; (5(5) $$\begin{array}{ccc}\hfill {\Phi }_C& =& {\gamma }_i\left(\theta \right){\Phi }_R+(1-{\gamma }_i\left(\theta \right)){\Phi }_{L}0.5\leqq {\gamma }_i\left(\theta \right)\leqq {\gamma }_i\left(0\right)\leqq 1.0\hfill \\ \hfill {\Phi }_C& =& {\gamma }_j\left(\theta \right)F_T+(1-{\gamma }_j\left(\theta \right)){\Phi }_{B}0.5\leqq {\gamma }_j\left(\theta \right)\leqq {\gamma }_j\left(0\right)\leqq 1.0\hfill \end{array}$$(5) ) continuous energy loss applied to charged-particle tracks between discrete interactions; (6(6) $$L{g}^{\left(n\right)}+R{g}^{\left(n\right)}=H{g}^{\left(n\right)}+Gg$$(6) ) total stopping power, which consists of soft Bremsstrahlung and collision loss terms; (7(7) $$f\times L{g}^{\left(n\right)}+f\times R{g}^{\left(n\right)}=f\times H{g}^{\left(n\right)}+Gg$$(7) ) collision loss determined by the restricted Bethe-Bloch stopping power with Sternheimer treatment of the density effect; (8(8) $${\Phi }_{\mathrm{unc}}\left(r\right)=exp\left(-{\sum }_t.|\mathit{r}-{\mathit{r}}_s|\right)/4\pi {|\mathit{r}-{\mathit{r}}_s|}^2$$(8) ) pair production; (9) Compton scattering; (10) coherent (Rayleigh) scattering, which can be included by means of an option; (11) the photoelectric effect.

Typical enhancements of EGS5 are described in the following sections.

3.2.1. Electron transport mechanics

A dual random hinge approach [13] has been adopted for the transport of electrons and positrons, in which energy loss and multiple elastic scattering are fully decoupled. The random multiple scattering hinge preserves near-second-order spatial moments of the transport equation over long step lengths. The hinge mechanics permits transport across boundaries between regions of differing media.

3.2.2. Binding effects in Compton interactions

Compton scattering is usually assumed that the electrons in the atomic clouds of the target atoms is free. Thus, it is ignored that the atomic electrons is bind to the nucleus. The effects, however, can have significant influence below tens of keV, particularly for the high-Z elements. Binding effects in the Compton model are included by taking into account the change in the total Compton cross section as well as the angular distribution of the emergent particles.

3.2.3. Doppler broadening and linearly polarized photon scattering

EGS5 can take into account the motion of the electrons in the atomic cloud. Because this motion constitutes a distribution of momenta, bound electrons ejected in Compton interactions emerge with a distribution of possible energies. This effect is usually called Doppler broadening. Polarized Compton and Rayleigh scattering are also included in EGS5.

Photon buildup factors and attenuation coefficients in water, iron, and lead were studied in the energy range 40–250 keV. The effects of binding, Doppler broadening, and polarization were pronounced in the lower energy range. As shown in Figure 4 on scattering photon spectrum from copper at the energy of 40 keV, the calculation with the effects reproduces the measurement almost perfectly.

Figure 4. Comparison among the calculation with and without the binding effect and the measurement on scattering photon spectrum from copper at the energy of 40 keV.

3.2.4. Electron impact ionization

K-shell electron-impact ionization (EII) is also included in EGS5. EII is treated as being a part of Møller scattering in EGS5, and neither the electron mean-free path nor the stopping power are modified by including EII.

Six different representations of the cross sections for EII are supplied: (1(1) $$\begin{array}{c}\hfill \frac{\mathit{dNND}\left(\mathit{r}\right)}{dr}=\frac{3}{2}·\frac{f_p}{1-f_p}·exp\left(-\frac{3}{2}·\frac{f_p}{1-f_p}·r\right),\end{array}$$(1) ) Casnati [37], (2(2) $$\begin{array}{ccc}& & \nabla ·\mathit{\Omega }\Phi (\mathit{r},\mathit{\Omega },E)+\sum _t(\mathit{r},E)\Phi (\mathit{r},\mathit{\Omega },E)\hfill \\ & & ={\int }_{E^{\text{'}}}d{E}^{\text{'}}{\int }_{\mathit{}{\mathbf{\Omega }}^{\text{'}}}d\mathit{}{\mathbf{\Omega }}^{\text{'}}\sum _S\mathit{r},E^{\text{'}}\to E,\mathit{}{\mathbf{\Omega }}^{\text{'}}\to \mathit{\Omega })\Phi (\mathit{r},\mathit{}{\mathbf{\Omega }}^{\text{'}},E^{\text{'}})\hfill \\ & & +S(\mathit{r},\mathit{\Omega },E)\hfill \end{array}$$(2) ) Kolbenstvedt (original) [38], (3(3) $$\begin{array}{ccc}& & {\Phi }_L\mu dy+{\Phi }_B\eta dx-{\Phi }_R\mu dy\hfill \\ & & -{\Phi }_T\eta dx-{\Sigma }_t{\Phi }_{C}dv=Q_{s}dv+Sdv\hfill \end{array}$$(3) ) Kolbenstvedt (revised) [39], (4(4) $$\begin{array}{c}{\Phi }_C=\alpha {\Phi }_R+(1-a){\Phi }_{L}0.5\leqq \alpha \leqq 1\hfill \\ {\Phi }_C=\beta {\Phi }_B+(1-\beta ){\Phi }_{T}0.5\leqq \beta \leqq 1\hfill \end{array}$$(4) ) Jakoby [40], (5(5) $$\begin{array}{ccc}\hfill {\Phi }_C& =& {\gamma }_i\left(\theta \right){\Phi }_R+(1-{\gamma }_i\left(\theta \right)){\Phi }_{L}0.5\leqq {\gamma }_i\left(\theta \right)\leqq {\gamma }_i\left(0\right)\leqq 1.0\hfill \\ \hfill {\Phi }_C& =& {\gamma }_j\left(\theta \right)F_T+(1-{\gamma }_j\left(\theta \right)){\Phi }_{B}0.5\leqq {\gamma }_j\left(\theta \right)\leqq {\gamma }_j\left(0\right)\leqq 1.0\hfill \end{array}$$(5) ) Gryzi´nski [41], and (6(6) $$L{g}^{\left(n\right)}+R{g}^{\left(n\right)}=H{g}^{\left(n\right)}+Gg$$(6) ) Gryzi´nski (relativistic) [42]. The user selects one of these six cross sections by specifying a PEGS5 input parameter.

After creation of a K-shell vacancy by EII, a K X-ray is emitted with the K-shell fluorescence yield, that is, the photoelectric effect, although the Auger electron or cascade particle is not treated in EII. For the case of emission, the X-ray energy is calculated from the ten K X-ray energies according to the relative yield of the K X-ray. The difference between the K-shell binding energy and the K X-ray energy is deposited locally.

3.3. PHITS [7]

PHITS (Particle and Heavy Ion Transport Code System) is a general-purpose Monte Carlo particle transport code written in FORTRAN. PHITS can deal with the transport of nearly all particles, including neutrons, protons, heavy ions, photons, and electrons, over wide energy ranges using several nuclear reaction models and nuclear data libraries. The geometrical configuration of the simulation can be set with general geometry (GG) or combinatorial geometry (CG). Various quantities, such as heat deposition, track length, and production yields, can be deduced from the PHITS simulation by using installed “tally” estimator functions. The code also has a function to draw two- and three-dimensional figures of the calculated results as well as the setup geometries, using the code ANGEL. PHITS can be executed on almost all computers (Windows, Mac, and Linux).

Figure 5 illustrates the current status of PHITS development and distribution. The main developments are carried out under the collaboration of JAEA, RIST, and KEK. JAEA is responsible for managing the entire project. In addition, there are five other contracts for PHITS development between JAEA and Chalmers University of Technology, Kyushu University, RIKEN, the French Atomic Energy Commission (CEA), and the Japanese Aerospace Exploration Agency (JAXA). The code has been distributed to many countries via RIST, the Data Bank of the Organization for Economic Co-operation and Development's Nuclear Energy Agency (OECE/NEA), and the Radiation Safety Information Computational Center (RSICC).

Figure 5. Current status of PHITS development and distribution.

Figure 6 summarizes the models used in PHITS for simulating nuclear and atomic collisions. The intra-nuclear cascade (INC) models, JAM [43] and INCL [44], are employed for simulating the dynamical stage of hadron-induced nuclear reactions for high- and intermediate-energy regions, respectively. Other INC models such as a modified version of Bertini and INC with emission of light fragment (INC-ELF) [45] can be also selected. For nucleus-induced reactions, a quantum molecular dynamics model JQMD [46] is employed because the consideration of nuclear forces between every combination of nucleons is crucial in the simulation. An evaporation and fission model GEM [47] is adopted for simulating the static stage. Nuclear and atomic data libraries are generally used for simulating low-energy neutron-induced nuclear reactions as well as photo- and electro-atomic interactions. These libraries were compiled on the basis of JENDL-4.0 and Lawrence Livermore Evaluated Electron Data Library (EEDL) [48].

Figure 6. List of models used in PHITS to simulate nuclear and atomic collisions.

PHITS has several important features, such as an event generator mode for low-energy neutron interaction, beam transport functions, a displacement per atom calculation function, and a microdosimetric tally function. Owing to these features, PHITS has been used for various applications such as accelerator design, medical and radiological protection, and cosmic-ray research. Some of these applications are introduced herein.

PHITS was extensively used in the shielding and neutron beam line design of the J-PARC project. Beam transport functions for simulating particle trajectories in magnet fields, gravity, T0-choppers, and neutron super mirrors are very useful in the design [49]. The calculations of effective doses and effective dose equivalents were carried out with PHITS by using the reference adult voxel phantoms. The event-generator mode was indispensable for this calculation. The calculated results were employed in the evaluation of the dose conversion coefficients by ICRP. The microdosimetric function of PHITS was used in the development of a new computational model for calculating the relative biological effectiveness (RBE)-weighted dose for charged particle therapy. The PHITS-based analytical radiation model in the atmosphere, PARMA [50], was established on the basis of PHITS simulation results for the atmospheric propagation of cosmic-rays. This model is currently employed for regulatory purposes to keep limit of the aircrew dose in Japan.

PHITS2 was officially released to the public in 2010. More than 800 persons in total from many different countries have been registered as PHITS2 users within three years, and the number is still rapidly increasing, by approximately 200 users each year. Thus, PHITS2 is one of the most successful nuclear-related codes developed mainly in Japan. A project for incorporating EGS5 into PHITS2 is currently in progress. The new PHITS2 has a potential to be an ultimate tool for simulating transport phenomena of all particles.

4. Discrete ordinate Sn transport method

The Sn transport method is a widely used discrete ordinate method to solve the Boltzmann transport equation. In this method, the spatial distribution of a neutron and/or photon flux is calculated by using the finite difference approximation, in which the calculation geometry is divided by spatial meshes into many small cells (or mesh intervals).

The whole energy region to calculate the neutron/photon flux is also divided into a number of energy groups, and the energy dependence of the neutron/photon reaction cross section is treated with a multi-group formulation. The Legendre expansion approximation is used to describe differentially scattering cross sections.

Anisotropy of the neutron/photon transport direction is described by a set of angular flux data given for a discrete direction. The name “Sn method” was derived from the angle-dependent flux description based on angular segmentation, in which the whole solid angle 4π is divided into N segments. An example of Sn angular segmentation is illustrated in Figure 7.

Figure 7. Example of Sn angular quadrature set: ω ≡ ω(μ,η). ω; direction vector of particle travel. μ η, ζ = (1-μ²-η²)^1/2; direction cosine of vector ω.

The accuracy of Sn calculation is determined by the adopted computational parameters: the number of energy groups, number of spatial meshes, number of angular direction points, and order of Legendre expansion.

The advantage of the Sn transport method is that the spatial distribution of energy-angle dependent flux can be calculated across the entire system.

However, adequate choices for the above-mentioned computational parameters are needed to achieve the required accuracy of flux solutions. For instance, rough spatial meshing sometimes causes negative flux results. Another problem characteristic of the Sn method is the so-called “ray effect”, i.e., calculation errors observed in large non-scattering regions, in which the radiation flux is enhanced along the direction of the Sn angular points.

In this section, Sn transport codes ANISN [51], DORT [52], and TORT [53] as well as the numerical methods applied in these codes are introduced. A computational technique developed to eliminate the ray effect will also be explained.

4.1. Progress in Sn transport codes

Discrete ordinate Sn method was originally developed by Bengt Carlson at Los Alamos National Laboratory (LANL) in 1953 [54]. Based on the theory of Carlson, the one-dimensional (1D) Sn transport code ANISN was developed at Oak Ridge National Laboratories (ORNL) in 1957. The ANISN code can be applied to a 1D slab/cylinder/sphere model.

In 1967, the two-dimensional (2D) Sn code DOT [55] was developed at ORNL, and the Sn transport method became an actual calculation method for practical shielding design. The final version of the DOT code series is the DORT code, which can treat discontinuous and triangular meshes. The DORT code can be applied to 2D XY, RZ, RΘ, and triangular geometry.

As for 3D Sn codes, the LANL THREETRAN [56], Japan's EMSAMBLE [57], and the ORNL TORT codes were developed in the 1980s. After undergoing various improvements, the TORT code became available for practical use in radiation transport calculation. Furthermore, recent progress in computers has largely expanded the practicable problem size of shielding calculation. However, use of the TORT code remains limited because of the enormous memory size and computation time required for an actual shielding design with 3D geometry. In the case of large-scale problems such as reactor in-vessel radiation calculation, the whole geometry is divided into several regions and radiation calculations for each region are performed separately.

The coupling of calculations in two different regions is performed by way of the so-called bootstrap method, in which the interior boundary flux picked up from the up-stream region is used as the boundary source of calculations for the down-stream region. The bootstrap method is applied to either 2D–3D or 3D–3D coupling. The auxiliary codes TORSED [58], TORSET [59], and VISA [58] are used for the procedure of bootstrapping.

The latest versions of the ANISN, DORT, and TORT codes and related auxiliary codes are included in the DOORS3.2a code system [60]. The calculation geometries addressed in Sn transport codes are illustrated in Figure 8.

Figure 8. Calculation geometries for one-, two-, and three-dimensional discrete ordinate transport codes.

4.2. Numerical method to solve Boltzmann equation

The steady-state Boltzmann transport equation is written as(2) $$\begin{array}{ccc}& & \nabla ·\mathit{\Omega }\Phi (\mathit{r},\mathit{\Omega },E)+\sum _t(\mathit{r},E)\Phi (\mathit{r},\mathit{\Omega },E)\hfill \\ & & ={\int }_{E^{\text{'}}}d{E}^{\text{'}}{\int }_{\mathit{}{\mathbf{\Omega }}^{\text{'}}}d\mathit{}{\mathbf{\Omega }}^{\text{'}}\sum _S\mathit{r},E^{\text{'}}\to E,\mathit{}{\mathbf{\Omega }}^{\text{'}}\to \mathit{\Omega })\Phi (\mathit{r},\mathit{}{\mathbf{\Omega }}^{\text{'}},E^{\text{'}})\hfill \\ & & +S(\mathit{r},\mathit{\Omega },E)\hfill \end{array}$$(2)

where Φ(r,Ω,E) is the angular flux for direction Ω and energy E at position r, and Σt is the total (absorption plus total scattering) cross section. ∇·ΩΦ indicates the divergence of the angular flux; Σt(r,E)Φ  denotes the decrease in flux Φ(r,Ω,E) due to absorption or scattering. Σs is the differential scattering cross section referring to the initial energy E' and direction Ω', as well as the final energy E and direction Ω after scattering. The term ∫∫ΣsΦ indicates the scattering source, which corresponds to the increase in flux Φ(r,Ω,E) due to down-scattering from higher energy groups and up-scattering from lower energy groups. In the form of a multi-group formulation, the term includes a contribution of within-group scattering. S(r,Ω,E) is a fixed source given as a calculation condition of the problem to be solved.

Figure 9 illustrates simplified two-dimensional model of the finite difference approximation method adopted in the Sn transport codes. The neutron/photon balance for an individual energy group and direction in an individual cell with a minute volume dv = dxdy is written as (3) $$\begin{array}{ccc}& & {\Phi }_L\mu dy+{\Phi }_B\eta dx-{\Phi }_R\mu dy\hfill \\ & & -{\Phi }_T\eta dx-{\Sigma }_t{\Phi }_{C}dv=Q_{s}dv+Sdv\hfill \end{array}$$(3) where the first four terms on the left indicate the net leakage, which is equivalent to the cell-integrated divergence of flux. The term –ΣtΦCdv indicates the decrease due to nuclear or electromagnetic interaction. The term Qsdv indicates the increase due to the scattering source. Qs is calculated by performing the integration ∫∫ΣsΦ in Equation (2)(2) $$\begin{array}{ccc}& & \nabla ·\mathit{\Omega }\Phi (\mathit{r},\mathit{\Omega },E)+\sum _t(\mathit{r},E)\Phi (\mathit{r},\mathit{\Omega },E)\hfill \\ & & ={\int }_{E^{\text{'}}}d{E}^{\text{'}}{\int }_{\mathit{}{\mathbf{\Omega }}^{\text{'}}}d\mathit{}{\mathbf{\Omega }}^{\text{'}}\sum _S\mathit{r},E^{\text{'}}\to E,\mathit{}{\mathbf{\Omega }}^{\text{'}}\to \mathit{\Omega })\Phi (\mathit{r},\mathit{}{\mathbf{\Omega }}^{\text{'}},E^{\text{'}})\hfill \\ & & +S(\mathit{r},\mathit{\Omega },E)\hfill \end{array}$$(2) , including revival by within-group scattering. The sum of Qs and the fixed source Sdv is the internal generation of neutrons/photons.

Figure 9. Particle (neutron/photon) balance.

Note: ${\Phi }_L\mu dy+{\Phi }_B\eta dx-{\Phi }_R\mu dy-{\Phi }_T\eta dx-{\sum }_t{\Phi }_{C}dv=Q_{s}dv+Sdv$. Φ_L; left boundary flux; Φ_B; bottom boundary flux; Φ_R; right boundary flux; Φ_T; top boundary flux; μ,η; direction cosine; Φ_C; flux at the center; Σ_t; total cross section; Q_s; scattering source; S; fixed source.

In each iteration of the Sn transport calculation, boundary fluxes ΦR and ΦT as well as the cell average flux ΦC are obtained by using the newest values of boundary fluxes ΦL and ΦB, which are the results of upstream cell calculations.

To obtain the solution of three unknown quantities from a single equation, the flux is extrapolated by assuming adequate within-cell distribution, as shown in Figure 10. In XY 2D geometry cases, for example, flux extrapolation is performed by using parameters α and β as follows. (4) $$\begin{array}{c}{\Phi }_C=\alpha {\Phi }_R+(1-a){\Phi }_{L}0.5\leqq \alpha \leqq 1\hfill \\ {\Phi }_C=\beta {\Phi }_B+(1-\beta ){\Phi }_{T}0.5\leqq \beta \leqq 1\hfill \end{array}$$(4)

Figure 10. Various models for within-cell flux extrapolation along X-axis. (1)  α = 0.5, linear model; (2)  α  = 1.0, step model; (3) 0.5 ≦ α ≦ 1, weighted difference model; (4) 0.5 ≦ γ(θ) ≦ α,  θ weighted difference model.

A schematic diagram of the assumptions made for a 2D flux distribution is illustrated in Figure 11. The extrapolation models adopted in Sn transport codes are (1(1) $$\begin{array}{c}\hfill \frac{\mathit{dNND}\left(\mathit{r}\right)}{dr}=\frac{3}{2}·\frac{f_p}{1-f_p}·exp\left(-\frac{3}{2}·\frac{f_p}{1-f_p}·r\right),\end{array}$$(1) ) linear, (2(2) $$\begin{array}{ccc}& & \nabla ·\mathit{\Omega }\Phi (\mathit{r},\mathit{\Omega },E)+\sum _t(\mathit{r},E)\Phi (\mathit{r},\mathit{\Omega },E)\hfill \\ & & ={\int }_{E^{\text{'}}}d{E}^{\text{'}}{\int }_{\mathit{}{\mathbf{\Omega }}^{\text{'}}}d\mathit{}{\mathbf{\Omega }}^{\text{'}}\sum _S\mathit{r},E^{\text{'}}\to E,\mathit{}{\mathbf{\Omega }}^{\text{'}}\to \mathit{\Omega })\Phi (\mathit{r},\mathit{}{\mathbf{\Omega }}^{\text{'}},E^{\text{'}})\hfill \\ & & +S(\mathit{r},\mathit{\Omega },E)\hfill \end{array}$$(2) ) step, (3(3) $$\begin{array}{ccc}& & {\Phi }_L\mu dy+{\Phi }_B\eta dx-{\Phi }_R\mu dy\hfill \\ & & -{\Phi }_T\eta dx-{\Sigma }_t{\Phi }_{C}dv=Q_{s}dv+Sdv\hfill \end{array}$$(3) ) weighted difference, and (4(4) $$\begin{array}{c}{\Phi }_C=\alpha {\Phi }_R+(1-a){\Phi }_{L}0.5\leqq \alpha \leqq 1\hfill \\ {\Phi }_C=\beta {\Phi }_B+(1-\beta ){\Phi }_{T}0.5\leqq \beta \leqq 1\hfill \end{array}$$(4) ) θ weighted difference models.

Figure 11. Schematic illustration of weighted/θ-weighted difference model to determine Φ_C, Φ_R-, and Φ_T.

The linear model makes the most likely assumption, in which ΦC is the simple average of boundary fluxes at both sides of the cell. However, this model sometimes results in a negative flux. Although the step model effectively avoids negative results, it is seldom used because of its poor accuracy.

The weighted difference model is a more practical one, which provides non-negative results close to the linear model assumption by choosing adequate values for “α” and “β”. If the coarseness of the spatial meshing is inadequate, parameters “α” and “β” will be set to near 1.0. In this case, the result of the weighted difference model becomes closer to that of the step model.

The θ weighted difference model achieves results that more closely approximate linear results. In this model, adjustable spatial extrapolation parameters γ_i(θ) and γ_j(θ) are used instead of α and β (5) $$\begin{array}{ccc}\hfill {\Phi }_C& =& {\gamma }_i\left(\theta \right){\Phi }_R+(1-{\gamma }_i\left(\theta \right)){\Phi }_{L}0.5\leqq {\gamma }_i\left(\theta \right)\leqq {\gamma }_i\left(0\right)\leqq 1.0\hfill \\ \hfill {\Phi }_C& =& {\gamma }_j\left(\theta \right)F_T+(1-{\gamma }_j\left(\theta \right)){\Phi }_{B}0.5\leqq {\gamma }_j\left(\theta \right)\leqq {\gamma }_j\left(0\right)\leqq 1.0\hfill \end{array}$$(5)

θ (0 ≦ θ ≦ 1) is a user input parameter to adjust the shape of the within-cell flux distribution. When θ is set to 0.0, the θ weighted difference model becomes equivalent to the weighted difference model. In general cases, a larger θ near 1.0 is chosen so as to obtain results that more closely approximate linear (i.e., more accurate) results. The recommended value of θ for use in DORT/TORT codes is 0.9.

4.3. Acceleration of flux convergence

Because convergence of transport calculations based on simple flux iteration is usually very slow, acceleration of flux convergence is indispensable for Sn transport codes. In this section, a simple “Groupwise rebalancing method” is introduced to briefly explain the acceleration technique. Practical methods of acceleration adopted in Sn transport codes are also introduced.

4.3.1 Groupwise rebalancing method

In the n-th iteration of transport calculation for energy group “g”, the balance equation is written as(6) $$L{g}^{\left(n\right)}+R{g}^{\left(n\right)}=H{g}^{\left(n\right)}+Gg$$(6)

where

Lg: leakage from the geometry boundary

Rg: removal by scattering and absorption

Hg: internal source corresponding to revival by within-group scattering

Gg: external source (fixed source, inflow from geometry boundary, generation due to scattering from the other energy groups)

The value of the internal source Hg⁽ⁿ⁾, in spite of being required as an input for Equation (6)(6) $$L{g}^{\left(n\right)}+R{g}^{\left(n\right)}=H{g}^{\left(n\right)}+Gg$$(6) , cannot be determined before the process of the n-th iteration. In the actual calculation, the preceding iteration result Hg^(n–1) is used instead of Hg⁽ⁿ⁾; (7) $$L{g}^{\left(n\right)}+R{g}^{\left(n\right)}=H{g}^{(n-1)}+Gg$$(7)

The inconsistency due to the difference between Hg^(n–1) and Hg⁽ⁿ⁾ delays the flux convergence. In the groupwise rebalancing method, this inconsistency is corrected by adjusting the flux result Φg⁽ⁿ⁾  by an appropriate flux multiplication factor “f”, which is determined so that the modification of flux data (i.e., Φg⁽ⁿ⁾→f Φg⁽ⁿ⁾) satisfies the following balance equation. (7) $$f\times L{g}^{\left(n\right)}+f\times R{g}^{\left(n\right)}=f\times H{g}^{\left(n\right)}+Gg$$(7)

Flux acceleration is thus achieved by performing the above procedure in each iteration of the transport calculation.

4.3.2 Partial current rebalancing method

Partial current rebalancing acceleration is performed in the same manner as the groupwise rebalancing method. The procedure of determining the multiplication factor f is performed for each mesh interval so that the balance equation is satisfied in every mesh interval. This method is much more effective than the simple groupwise rebalancing method.

4.3.3 Diffusion synthetic acceleration [61]

The diffusion synthetic method is based on a procedure deferent from that of the groupwise rebalancing or partial current rebalancing method. The scattering source term, a necessary input for the mesh-wise balance equation, is prepared in parallel by performing a modified diffusion calculation in which either the diffusion coefficient, removal cross section, or external source is modified to obtain a precise source distribution that is consistent with the transport solution.

The diffusion synthetic method is as effective as the partial current rebalancing method. Figure 12 compares the convergence speeds obtained in 20 × 20 mesh XY geometry calculations with three different acceleration methods.

Figure 12. Convergence of flux solution observed in 20 × 20 mesh problem with a source region at the corner mesh interval.

4.3.4. Ray effect

When 2D or 3D Sn transport code is applied to a problem with large non-scattering regions, the flux solution is strongly affected by a computational error called the “ray effect”. A typical ray effect is observed in the analysis of a detector response apart from the source material and the analysis of radiation skyshine. In this section, several computational techniques to eliminate the ray effect are introduced.

1. Last collision source method

In the analysis of a detector response outside the geometry containing the source, the radiation transmission in the surrounding air/void region is calculated separately from the Sn transport for the source geometry.

To eliminate the ray effect, the air/void transmission is calculated on the basis of an analytical method. The widely used FALSTF [62] and SPACETRAN [63] codes are developed for calculating the air/void region.

The FALSTF code is based on the last collision source method. Prior to the calculation, the radiation flux and moment inside the geometry are obtained with the DORT code. The FALSTF code calculates the direct transportation of uncollided flux from mesh points inside the DORT calculation to the detector location. The detector response is calculated by integrating the uncollided flux contribution over the DORT calculation geometry.

The SPACETRAN code is based on a simple geometric method. Prior to SPACETRAN calculation, the energy- and angle-dependent radiation leakage from the geometry surface is prepared by DOT/DORT transport calculation. The radiation flux at the detector location is obtained by integrating the direct contribution of leakage radiation over the surface of the DOT/DORT calculation geometry.

(2)	First collision source method

The calculation geometry of skyshine analysis usually consists of a large air region and a small source region at the corner. When the radiation is calculated with the ordinary Sn transport procedure, in which a mesh-by-mesh sweep over angular fluxes is conducted from the source region to the geometry boundaries, the calculated result will be strongly affected by the ray effect. In the case of skyshine analyses, the first collision source method is effective in reducing the ray effect. In this method, the uncollided flux Φ_unc at an arbitrary point r in the calculation geometry is written as (8) $${\Phi }_{\mathrm{unc}}\left(r\right)=exp\left(-{\sum }_t.|\mathit{r}-{\mathit{r}}_s|\right)/4\pi {|\mathit{r}-{\mathit{r}}_s|}^2$$(8) where Σt is the total cross section of air, rs is the coordinate of the source, and |r – rs| is the distance from the source. The first collision source is the product of the uncollided flux and differential scattering cross section. Anisotropy of the first collision source is expressed in the form of a spherical harmonic expansion.

GRTUNCL [64] is an analytic first collision source code developed for DORT RZ geometry calculations. This code can treat a point source or a ring-shaped source as an actual source from which neutrons and/or photons are emitted. The first collision source corresponding to a disk source or cylindrical source can also be prepared by performing a series of GRTUNCL calculations with ring sources differing in size and location.

As codes for 3D first collision sources, GRTUNCL-3D [65] and FNSUNCL3 [66] were developed for TORT XYZ geometry calculations. Furthermore, GRTUNCL-3D has reportedly been applied to TORT RΘZ cylindrical geometry [67].

Extensive development of the first collision source method enabled the Sn transport method to treat skyshine analyses with various source geometry conditions.

5. Simplified calculation method

As mentioned in the previous sections, with the dramatic increase in computer power a detailed calculation technique such as the Monte Carlo method can now be applied to entire facilities. Simplified calculation methods, however, remain widely used in the various facilities to survey parameters before designing the final shielding and cross-checking the results by detailed calculation methods. Basic parameters such as gamma-ray buildup factors can be improved by using detailed calculation methods. However, these studies are mainly conducted in Japan. A special committee on Shielding Safety Evaluation Methods and Related Data in Nuclear Facilities was organized under the AESJ from 1995 to 2000. In this committee, various basic parameters such as gamma-ray buildup factors, response functions for skyshine, dose attenuation of neutrons, and databases for duct streaming were studied for simplified calculation methods. These activities were presented at the Nineth International Conference on Radiation Shielding (ICRS-9) held at Tsukuba in 2000 [68--76]. In this section, Japanese research on basic parameters such as gamma-ray buildup factors and response functions for skyshine as well as on simple calculation codes is introduced.

5.1. Basic parameters

5.1.1. Gamma-ray buildup factors

The gamma-ray buildup factor is defined as the ratio of the total value of a specified radiation quantity at point to the value from radiations reaching the point without undergoing a collision.

The current standard used across the world, ANSI/ANS-6.4.3 [6], was prepared by the working group ANS-6.4.3, which was formed in 1980 to develop a reliable and comprehensive data-set that reflected the experience of the accident at Three Mile Island. The codes used were ADJMOM-I [77](moments method), PALLAS [78,79](discrete ordinate-integral transport), and ASFIT [80] (discrete ordinate-integral transport). The data selected for the standard are mostly from moment method calculations. The materials are Be, B, C, N, O, Mg, Al, Si, P, Ar, K, Ca, Fe, Cu, water, air, and concrete. Data for the remaining materials, such as Mo, Sn, La, Gd, W, Pb, and U, were taken from PALLAS calculations. The geometric progression (G-P) form [81–83], developed by Harima et al., was selected as the best available form and included as part of ANSI/ANS-6.4.3. Contributions from Japan to this working group were explained in detail in elsewhere [84].

After the publication of ANSI/ANS-6.4.3, systematic studies were continued in Japan by the special committee organized in the AESJ. Following the committee, Shimizu et al. generated a new data-set of buildup factors by the invariant embedding (IE) method [85]. Specific features of the data-set as compared with the ANSI/ANS-6.4.3 standard, the data-set has the following unique features: (1(1) $$\begin{array}{c}\hfill \frac{\mathit{dNND}\left(\mathit{r}\right)}{dr}=\frac{3}{2}·\frac{f_p}{1-f_p}·exp\left(-\frac{3}{2}·\frac{f_p}{1-f_p}·r\right),\end{array}$$(1) ) data corresponding to the effective dose with AP and ISO irradiation geometries, as well as ambient dose equivalents are added; (2(2) $$\begin{array}{ccc}& & \nabla ·\mathit{\Omega }\Phi (\mathit{r},\mathit{\Omega },E)+\sum _t(\mathit{r},E)\Phi (\mathit{r},\mathit{\Omega },E)\hfill \\ & & ={\int }_{E^{\text{'}}}d{E}^{\text{'}}{\int }_{\mathit{}{\mathbf{\Omega }}^{\text{'}}}d\mathit{}{\mathbf{\Omega }}^{\text{'}}\sum _S\mathit{r},E^{\text{'}}\to E,\mathit{}{\mathbf{\Omega }}^{\text{'}}\to \mathit{\Omega })\Phi (\mathit{r},\mathit{}{\mathbf{\Omega }}^{\text{'}},E^{\text{'}})\hfill \\ & & +S(\mathit{r},\mathit{\Omega },E)\hfill \end{array}$$(2) ) the distance from the source is extended up to 100 mean free paths (mfp) from 40 mfp in the ANS data; (3(3) $$\begin{array}{ccc}& & {\Phi }_L\mu dy+{\Phi }_B\eta dx-{\Phi }_R\mu dy\hfill \\ & & -{\Phi }_T\eta dx-{\Sigma }_t{\Phi }_{C}dv=Q_{s}dv+Sdv\hfill \end{array}$$(3) ) the IE method and the photon cross section PHOTX are used consistently for all materials; and (4(4) $$\begin{array}{c}{\Phi }_C=\alpha {\Phi }_R+(1-a){\Phi }_{L}0.5\leqq \alpha \leqq 1\hfill \\ {\Phi }_C=\beta {\Phi }_B+(1-\beta ){\Phi }_{T}0.5\leqq \beta \leqq 1\hfill \end{array}$$(4) ) the treatment of Bremsstrahlung is improved by using the Bremsstrahlung production data calculated by the EGS4 code. As an example of the improvement, Table 2 [86] shows the maximum difference in the ratio of buildup factor with and without Bremsstrahlung between 0.5 and 40 mfp for materials not considering Bremsstrahlung in the ANSI/ANS-6.4.3. The data-set has been evaluated by the working group on gamma-ray buildup factors of the subcommittee on radiation shielding, and prepared as standard data in the AESJ. This activity was reported at the 12th International Conference on Radiation Shielding held in Nara in 2012 [87].

Table 2. Effect of Bremsstrahlung to gamma-ray buildup factors. (Maximum ratio of buildup with and without Bremsstrahlung between 0.5 and 40 mfp for each energy and material)

Material/Energy	Be	B	C	N	O	Na	Mg	Al	Si
15.0 MeV	1.072	1.082	1.094	1.098	1.113	1.186	1.213	1.242	1.268
10.0 MeV	1.043	1.046	1.051	1.052	1.057	1.087	1.097	1.108	1.118
8.0 MeV	1.034	1.035	1.037	1.037	1.04	1.059	1.065	1.071	1.077
6.0 MeV	1.024	1.024	1.024	1.024	1.025	1.035	1.039	1.042	1.045
5.0 MeV	1.02	1.019	1.019	1.019	1.019	1.026	1.028	1.03	1.032
4.0 MeV	1.016	1.014	1.014	1.013	1.013	1.017	1.019	1.02	1.021
3.0 MeV	1.012	1.01	1.009	1.009	1.008	1.01	1.011	1.012	1.012
2.0 MeV	1.008	1.006	1.005	1.005	1.004	1.005	1.005	1.005	1.006
1.5 MeV	1.005	1.004	1.003	1.003	1.003	1.003	1.003	1.003	1.003
Material/Energy	P	S	Ar	K	Fe	Cu	Water	Concrete	Air
15.0 MeV	1.297	1.327	1.378	1.433	1.808	2.044	1.093	1.203	1.098
10.0 MeV	1.127	1.138	1.152	1.169	1.276	1.328	1.047	1.094	1.049
8.0 MeV	1.082	1.088	1.095	1.106	1.163	1.188	1.032	1.063	1.034
6.0 MeV	1.048	1.05	1.053	1.058	1.086	1.095	1.02	1.038	1.021
5.0 MeV	1.033	1.035	1.037	1.04	1.057	1.063	1.015	1.027	1.016
4.0 MeV	1.022	1.023	1.024	1.026	1.036	1.038	1.01	1.018	1.011
3.0 MeV	1.013	1.013	1.014	1.014	1.019	1.02	1.006	1.011	1.007
2.0 MeV	1.006	1.006	1.006	1.006	1.008	1.008	1.003	1.005	1.003
1.5 MeV	1.003	1.003	1.003	1.003	1.004	1.004	1.002	1.003	1.002

5.1.2. Response functions for skyshine

The line- or conical-beam response function (LBRF/CBRF) is the air-scattered dose from radiation emitted at a fixed angle (LBRF) or fixed opening angle (CBRF) plotted as a function of the distance separating the source and detector. LBRF was first presented by Lynch et al. [88] as air-scattered flux for gamma rays with the Monte Carlo method.

Gamma-ray LBRF and CBRF were calculated by using the EGS4 Monte Carlo code and studied as part of the activities of the special committee. Calculated databases of LBRF and CBRF were summarized and published as the KEK Report, together with their four-parameter empirical formulas [89].

Neutron CBRF up to 3 GeV was also calculated and studied as a part of the activities of the special committee. The final version of CBRF was calculated by using PHITS and used for the basic parameters of the SHINE3 code [90].

5.2. Simple calculation code

Despite the progress of research on basic parameters in Japan, simple calculation codes, especially the most widely used point kernel codes for gamma-rays, are based on those developed in USA such as QAD [11], G33 [91], and SKYSHINE [92]. Attempts have been made to update the basic parameters obtained from research in Japan into such codes especially by using GP parameters such as QAD-CGGP [93]. In this case, the right of improvement belongs to the original authors of the codes. Ideally, creating and maintaining point kernel codes for gamma-rays in Japan to utilize the active research on basic parameters is desirable.

In contrast, simple calculation codes for accelerators have been developed and maintained in Japan at KEK and JAEA. BULK-I [94] is the bulk shielding calculation code for 50-500MeV proton accelerator facilities, which are mainly installed in cancer treatment facilities that use protons. DUCT-III [95] is the duct streaming code for gamma-rays below 10 MeV and neutrons below 3GeV. As mentioned previously, SHIN3 is the skyshine calculation code for neutrons up to 3GeV.

6. Shielding experiments

Many types of shielding experiments have been carried out on radioisotope irradiation facilities, fission reactors, fusion research fields, and accelerator facilities. The experimental data have been used for various purposes: to deduce and validate empirical formulas, to evaluate parameters for analytical equations, and to validate shielding design codes. Most empirical formulas on streaming phenomena and their parameters, i.e. albedo data, have been deduced from the results of streaming experiments. Buildup factors and attenuation lengths used in the analytical equations have also been evaluated by experimental results. Experimental data have been used to validate the analytical equations for skyshine phenomena as well. Recently, most shielding experiments have been used for benchmarking on shielding design codes with the effect of dramatically improving the codes.

Benchmark experiments are categorized as follows, and their scale and required precision vary depending on the category.

1. Integral experiments: Experiments using a simple experimental setup, with a simple material and composition for precise evaluation of cross sections.

2. Basic benchmark experiments: Experiments to validate codes, empirical formulas, and analytical equations.

3. Engineering benchmark experiments: Experiments with typical geometrical structures to verify the reliability of overall shielding designs.

4. Mock-up experiments: Experiments to verify the accuracy and adequacy of engineering designs.

In benchmark experiments, all experimental conditions and data required for analyses by the codes must be clear and available. Source, geometrical, and boundary conditions are the most important data for benchmarking. The quality of the benchmark experiment depends on the physical quantity of the measured data. In many experiments, data such as dose and reaction rate as a function of position are generally measured by dosimeters and detectors with ease. More useful and effective for benchmarking the code are data such as energy spectra as a function of position and energy, which are measured in integral and basic benchmark experiments. Additionally, angle-dependent energy spectra, which are occasionally measured in integral experiments, are the most useful and effective for benchmarking.

Benchmark experiments on fission reactors and radioisotope irradiation facilities were compiled and released twice by the shielding committee of AESJ as shielding benchmark problems to validate the codes, models, and related nuclear data systematically [96,97]. The benchmark experiments have also been compiled as SINBAD (the radiation shielding experiments database) since 1980, with the newest version released in 2011 [98]. SINBAD is maintained as a basis for development and validation of computer code, model, and nuclear data. The collection of experimental data is jointly managed by Organization for Economic Co-operation and Development, Nuclear Energy Agency, Data Bank (OECD/NEA/DB), and RSICC. The number of data currently contained in SINBAD are 45, 31, and 23 experiments on reactor shielding, fusion neutronics, and accelerator shielding, respectively. The compiled data are found in

1. SINBAD reactor shielding benchmark experiments

2. SINBAD fusion neutronics benchmark experiments

3. SINBAD accelerator shielding benchmark experiments

Of the reactor shielding benchmark experiments, reactor dosimetry experiments carried out at the H.B. Robinson II and VENUS-3 facilities in USA and the Balakovo Unit 3 facility in Russia are reviewed and covered as excellent benchmark problems. The following experiments are also reviewed: deep penetration experiments at the Tower Shielding Facility of the ORNL, at the ASPIS shielding facility of the Winfrith site, and at the Yayoi of the University of Tokyo; and the streaming experiment at IRI of the Technical University of Budapest. The Japanese–American Shielding Program of Experimental Research (JASPER) program, under the collaboration of Japan and USA, has provided excellent shielding data for benchmarking on liquid metal fast breeder reactors.

Many of the integral experiments reviewed are infusion neutronics benchmark experiments, i.e., tests carried out at OKTAVIAN of the Osaka University, Fusion Neutronics Source (FNS) of JAEA, Generatore di Neutroni Frascati (FNG) of the Italian National Agency for New Technologies, Energy and Sustainable Economic Development (ENEA), and Technische Universităt Dresden (TUD) of Germany. Deep penetration, streaming, and skyshine experiments carried out at these facilities as basic benchmark experiments are also reviewed. Recently, mock-up experiments are also being conducted by the facilities under the collaborative framework of the ITER project [99,100], although the review has not yet been finalized for SINBAD.

Of the accelerator shielding benchmark experiments, source term and neutron transmission experiments are mainly covered and reviewed. Neutron energy spectra from various thick targets bombarded by various charged particles were measured at the National Superconducting Cyclotron Laboratory facility (NSCL ) of the Michigan State University (MSU) and the Heavy Ion Medical Accelerator in Chiba (HIMAC) of the National Institute of Radiological Sciences (NIRS). Transmission experiments were carried out at ISIS pulsed neutron and muon source of the Rutherford Appleton Laboratory and the Takasaki Ion Accelerators for Advanced Radiation Application (TIARA) of JAEA. Various types of experiment carried out at the CERN-EU High Energy Reference Field Facility (CERF) are also covered in SINBAD. Recently, some shielding experiments have been carried out and reported at the Brookhaven National Laboratory, the Fermi National Accelerator Laboratory, and Research Center for Nuclear Physics of the Osaka University (RCNP) [101–103].

Review activity of data in SINBAD is ongoing for future inclusion in the revised SINBAD. The review activity includes the completeness and consistency of information on experiments and evaluates the experimental sources of uncertainty.

Benchmark experiment databases such as SINBAD will become increasingly important for the development and maintenance of shielding codes, models, and nuclear data. Because code verification is carried out on the basis of benchmark experiment data, the benchmark experiments is rigorously reviewed under strict criteria, and, if possible, under international activity. An international committee is established to authorize the review criteria on experiments. The database needs to cover a wide energy region and various shielding items such as streaming and skyshine for reactor, fusion, and accelerator applications, as well as all categories of benchmark experiments. It is better that more data are collected to satisfy the database requirements. Some types of benchmark experiment may be proposed by a scientific authority such as an international committee.

6. Conclusion

In this review, the “Technological Roadmap in the Radiation Shielding Field,” prepared by the “Radiation Shielding” research committee established in the AESJ, is first presented for considering future prospects, technological levels that need to be confirmed in light of current needs, and remaining issues that need to be resolved for technological development. Subsequently, the progress and prospects of individual calculation methods for radiation shielding, such as the Monte Carlo, discrete ordinate Sn transport, and simplified calculation methods, are presented. Finally, shielding experiments used to establish benchmarks for shielding calculation methods are presented.

Activities of the standardization of shielding codes recently appear for licensing in Japan. Currently, old shielding design codes developed and imported a few decades ago have been used mainly for licensing in Japan, although many codes have been modified and newly developed. Thus, the newest knowledge has not been applied for licensing except for some special cases. Moreover, code copyright and export control problems have come to the forefront in the last few years. To overcome these problems, it is better that a protocol to accept new codes is defined and stated clearly by an authority. The protocol may include processes to develop codes, evidence of reliability of the codes provided by the developers, an endorsement of the codes by authorities, and a demonstration that the code is standardized. Additionally to guarantee reliability, validation and verification of the codes may be carried out in the processes of development and maintenance.

Acknowledgements

The authors wish to thank Dr Y. Nagaya and T. Sato of JAEA and Prof. Y. Namito of KEK for their supports of description on the Monte Carlo methods.