Neutronic calculations of the China dual-functional lithium-lead test blanket module with the parallel discrete ordinates code Hydra

Guang-Chun Zhang· Jie Liu · Liang-Zhi Cao· Hong-Chun Wu ·Xian-Bao Yuan

Abstract The China dual-functional lithium-lead test blanket module (DFLL-TBM) is a liquid LiPb blanket concept developed by the Institute of Nuclear Energy Safety Technology of the Chinese Academy of Sciences for testing in ITER to validate relevant tritium breeding and shielding technologies. In this study, neutronic calculations of DFLL-TBM were carried out using a massively parallel three-dimensional transport code, Hydra, with the Fusion Evaluated Nuclear Data Library/MG. Hydra was developed by the Nuclear Engineering Computational Physics Lab based on the discrete ordinates method and has been devoted to neutronic analysis and shielding evaluation for nuclear facilities. An in-house Monte Carlo code(MCX) was employed to verify the discretized calculation model used by Hydra for the DFLL-TBM calculations.The results showed two key aspects:(1)In most material zones,Hydra solutions are in good agreement with the reference MCX results within 1%, and the maximal relative difference of the neutron flux is merely 3%, demonstrating the correctness of the calculation model; (2) while the current DFLL-TBM design meets the operation shielding requirement of ITER for 4 years, it does not satisfy the tritium self-sufficiency requirement. Compared to the two-step approach, Hydra produces higher accuracies as it does not rely on the homogenization technique during the calculation process. The parallel efficiency tests of Hydra using the DFLL-TBM model also showed that this code maintains a high parallel efficiency on O(100)processors and,as a result, is able to significantly improve computing performance through parallelization. Parameter studies have been carried out by varying the thickness of the beryllium armor layer and the tritium breeding zone to understand the influence of the beryllium layer and breeding zone thickness on tritium breeding performance. This establishes a foundation for further improvement in the tritium production performance of DFLL-TBM.

Keywords Discrete ordinates method · DFLL-TBM ·Neutronic analysis · Tritium breeding performance

1 Introduction

The China dual-functional lithium-lead test blanket module (DFLL-TBM) [1, 2] is a liquid LiPb blanket concept for ITER developed by the Institute of Nuclear Energy Safety Technology (INEST) of the Chinese Academy of Sciences. It was designed to demonstrate the integrated technologies of the helium single coolant (SLL) blanket and the He/LiPb dual coolant (DLL) blanket. To validate the feasibility of DFLL-TBM, a series of theoretical analyses and numerical experiments [3-5] have recently been carried out to assess its thermal and mechanical properties.

Neutronics analyses of the TBM play a vital role in the design, construction, and experimentation of ITER. A previous study [6] proposed a two-step calculation scheme combining the method of characteristic(MOC)and the simplified spherical harmonics (SPN) method to carry out neutronic calculations for the DFLL-TBM.In this twostep approach, the three-dimensional (3D) TBM model is axially divided into several planes, and planes with unique geometry and material distributions are selected as the‘‘typical planes’’ on which the MOC calculation is conducted to determine the detailed fine-group flux distributions. The homogenization technique is subsequently utilized to condense the fine-group cross sections.Lastly,a 3D SPNcalculation is undertaken to obtain the coarsegroup flux profile throughout the TBM model. The twostep scheme is a compromise between solution accuracy and calculation efficiency. Therefore, this two-step approach yields results with relatively large errors; this is the underlying rationale driving the pursuit of a direct 3D deterministic neutronic solution for TBM problems.

In this work, particular attention is given to the nuclear assessment of DFLL-TBM, including tritium breeding and nuclear shielding performances, using a 3D deterministic transport method. Deterministic transport methods may be classified into two categories based on the approximations for discretizing the angular variables:the discrete ordinates(SN) method [7] and the spherical harmonics (PN) method[8, 9]. The PNmethod usually leads to sophisticated algebraic systems and is rarely used for practical 3D calculations.On the other hand,the discrete ordinates(SN)method is considered one of the most important and effective tools in the neutronic analyses field due to its computational appeal.The very first development of the SNmethod dates back to the 1950s;Carlson[7]described a finite difference approximation to the Boltzmann transport equation in which the angular variable is represented by piecewise continuous linear functions. Over recent decades, there have been tremendous research efforts to enhance solution accuracy [10] and iterative stability [11] and develop effective acceleration schemes[12]for the SNmethod.The effectiveness, efficiency and generality of the SNmethod rather quickly led to the development of several one-dimensional (1D) [13], two-dimensional (2D) [14], and 3D[15] SNcodes, some of which were heavily used in the nuclear analyses field [16, 17]. With rapidly evolving computational resources, there has been an increasing interest in developing massively parallel 3D SNcodes[18-20] and applying them to neutronic calculations for TBMs.Various parallel algorithms for the SNmethod have been proposed over the last decade. Among them, the Koch-Baker-Alcouffe (KBA) algorithm [21] is noted for its high efficiency for structured 3D mesh grids and has been used extensively in many modern parallel SNtransport codes.

Recently,a massively parallel 3D transport code,Hydra[22], based on the SNmethod has been developed by the Nuclear Engineering Computational Physics Lab (NECP)of Xi’an Jiaotong University. Hydra is characterized by a mixed programming of C++ and Fortran, domain decomposed parallel computation based on the KBA algorithm, and the utilization of various spatial discretization schemes. Hydra can also conduct criticality analyses for nuclear reactors, in addition to shielding analyses for nuclear facilities. Hydra has been verified by a variety of 3D benchmark problems[22,23],and the numerical results have demonstrated that this code possesses very high solution accuracies. Hydra substantially reduces the computational time by parallelizing the transport sweeping procedure,when compared to the well-known SNtransport code, TORT [15], developed by the Oak Ridge National Laboratory (ORNL). Currently, both the Cartesian and cylindrical meshes may be used by Hydra to discretize 3D models. In contrast, in the Denovo [24] code, only the Cartesian mesh discretization is implemented. As a comparison, the ATTILA [25] code is able to utilize unstructured tetrahedra meshes by employing linear discontinuous finite element method to discretize the spatial variables.

In this study, we applied the Hydra code to neutronic assessments for DFLL-TBM. The discretized calculation model used by Hydra was verified against an in-house Monte Carlo code, MCX, developed by the NECP laboratory. Various neutronic quantities were obtained including the neutron flux profile along the radial direction,tritium breeding ratio (TBR) and tritium production rate,displacement per atom (DPA), and hydrogen and helium production rates. We also compared the results of the onestep approach based on Hydra with those of the two-step method.The comparison results demonstrated that the onestep method yields solutions with higher accuracy as it avoids the homogenization process.A parameter study was carried out by varying the thickness of the tritium breeding zone, and the first wall beryllium layer to evaluate how tritium breeding performance was influenced by the beryllium layer and breeding zone thickness.This provides the foundation for further improvement of the performance of DFLL-TBM.

The remainder of this paper is organized into five sections. In Sect. 2, the methodologies of the 3D transport solver of Hydra used in this study’s simulations are briefly reviewed.Section 3 describes the geometry and calculation model of DFLL-TBM, and Sect. 4 details the results.Section 5 presents the parameter study results with respect to the tritium breeding zone and the Be armor thickness.The discussion and conclusions are provided in Sect. 6.

2 Methodologies

where φgmis the angular flux in direction Ωm=（μm，ηm，ξm）; Qgmis the group source including the scattering and external fixed sources; and M is the number of discretized angles.

Hydra attempts to discretize the system of equations of Eq. (2) with the difference method. The resultant discretized transport problem can be written in the operator notation as:

where Lg=（Ω·∇+Σg） is the transport operator; φgand Qgare the vectors of angular fluxes and group sources,respectively, for every angle in each spatial mesh.

The source iteration (SI) method, also known as the Richardson iteration scheme, is used in Hydra to solve the linear system of Eq. (3). In the SI method, the transport sweep process is carried out repeatedly to converge the group fluxes. In problems where the upscatter is present,iterations between thermal groups should be performed to converge upscattering group sources. The overall calculation procedure of the 3D transport solver of Hydra for fixed-source problems is illustrated in Fig. 1.

Currently, Hydra supports two different types of structured mesh configurations: the Cartesian and cylindrical mesh geometries as illustrated in Fig. 2. Based on the specific geometry of DFLL-TBM,only the Cartesian meshgrid was used to discretize the calculation model.Both the volumetric and surfaces sources can be set in the Hydra code to address different types of fixed-source problems.In this study, the surface source was used to simulate highenergy fusion neutrons emitted from the plasma zone.

The transport solution options of Hydra used in this work involve several important implementation details.More specific,the level-symmetric quadrature set,diamond difference (DD) and theta-weighted diamond (TWD) difference schemes,2D decomposition of 3D meshes,and the KBA parallel transport sweep algorithm were used excessively in neutronic calculations for DFLL-TBM and are discussed in Sects. 2.1-2.4.

2.1 Level-symmetric quadrature set

The SNmethod relies on the quadrature set to evaluate integrals over directional space. The quadrature set is a collection of discrete directions and associated weights.The angular flux is only evaluated in these directions, and the scalar flux is estimated by the weighted sums of the angular flux. Over the years, a variety of quadrature set configurations have been developed by imposing distinct constraints, such as achieving rational symmetry [15] or higher-order precision [26, 27]. Among them, the levelsymmetric quadrature set has been adopted and implemented in Hydra. The advantage of the level-symmetric quadrature set is that it is invariant with respect to a 90°axis rotation. The level-symmetric LS6quadrature set configuration of an octant is illustrated in Fig. 3.

2.2 Spatial discretization schemes

The diamond difference scheme is one of the simplest forms of spatial discretization. It approximates the angular flux over a Cartesian cell using a linear function between adjacent boundaries.This means the average angular flux is simply the linear average of the boundary values. The outgoing angular flux is then determined by extrapolating the incoming and average angular fluxes. Using the mesh shown in Fig. 4 as an example, the average angular flux φi，j，k，m，gfor the mesh cell （i，j，k） along direction, Ωm=（μm，ηm，ξm） in group, g, can be expressed as:

It should be noted that due to the extrapolation technique used to determine the outgoing flux shown in Eq. (5),a negative solution may arise if the optical thickness of the spatial cell is too large and consequently thwarts the iteration convergence.This problem may be remedied with the TWD scheme which uses the incoming fluxes to calculate weighting factors to permit positivity of the cell-centered and outgoing angular fluxes.The weighting factors a,b and c are calculated using the following system of equations:

where

The cell-centered and outgoing angular fluxes can then be calculated as:

Note that Hydra uses the same default value of θ=0.9,as used in the Denovo [24] and TORT [15] codes. The derivation of the spatial discretization schemes of Hydra is described in detail in Ref [22].

2.3 Spatial domain decomposition

The spatial domain decomposition (SDD) method [28]is a widely used approach to conduct distributed-memory parallel computing. The SDD’s typical procedure for a neutron transport problem involves subdividing the spatial domain of the problem into subdomains, conducting transport sweeps independently on each subdomain and communicating angular flux on subdomain boundaries whenever necessary. In general, 3D SDD is employed where the domain decomposition is performed not only in the radial direction, but also in the axial direction. However,the KBA parallel transport sweep algorithm requires a 2D decomposition of the 3D mesh, as illustrated in Fig. 5.This figure shows that the mesh decomposition is only conducted on the 2D X-Y plane.All meshes stacked in the axial direction reside in the same processor.

2.4 KBA parallel transport sweep algorithm

The main computational burden of a 3D transport code originates from repeated executions of the transport sweep process used to invert the transport operator of the neutron transport equation. Therefore, the transport sweep process must be parallelized to achieve a good computational performance. The KBA algorithm is a widespread acceptance technique for parallelizing the transport sweep process on orthogonal structured grids. It differs from other parallel computing algorithms such as the parallel block Jacobi (PBJ) method in that KBA is a direct inversion of the transport operator, while PBJ inverts the transport operator in an iterative manner.

Figure 6 illustrates an example of the KBA parallel transport sweep algorithm for the 3D problem in Fig. 5.Note that the blue cube denotes the smallest parallel block comprised of multiple computational mesh cells. The transport sweep process can be conducted simultaneously among certain computational blocks such as the three blocks shown in Fig. 6b and the six blocks shown in Fig. 6c. However,inside each parallel block the transport sweep is still executed serially. It can be seen from Fig. 6 that all parallel computational blocks form a 2D plane propagating from a corner of the 3D grids along the neutron streaming direction.For this reason,the KBA algorithm has also been referred as the wave front method in some references[29,30].

3 The DFLL-TBM geometry and calculation model

Figure 7 shows that the DFLL-TBM model consists of the first wall(FW),radial-poloidal stiffening plates(rpSPs)and toroidal-poloidal stiffening plates (tpSPs), LiPbchannels, and back planes (BPs). The FW is designed to withstand the heat flux from the plasma chamber and is covered by a 2-mm-thick beryllium layer.It is composed of three layers(FW1,FW2 and FW3),and the layout of these three layers is illustrated in Fig. 5.FW1 and FW3 are both made of the China low activation martensitic steel(CLAM)[31]. FW2 is a mixture of CLAM with 20.4 vol% and helium with 79.6 vol%.The BPs are designed as the strong back support for TBM and also act as a helium manifold.The BPs are made of two thick plates (BP1, BP4: 20 mm thick) serving as structure functions, and two intermediate thin plates (BP2, BP3: 10 mm thick) serving as flow separators.

The DFLL-TBM relies on the use of Pb-Li liquid eutectic alloy, both as tritium breeder and neutron multiplier. The rpSP and tpSP, shown in Fig. 7, separate the tritium breeder zone into six channels. To facilitate comparisons in the following calculations, the six LiPb channels were grouped into three zones, namely the LiPb2,LiPb2, and LiPb3 zones. The neutron flux and tritium productions were evaluated based on these three zone configurations.

Table 1 lists the materials and dimensions for the DFLL-TBM[1].Detailed descriptions on the geometry and material configuration are provided in Refs. [1, 2].

The discretized calculation model of the DFLL-TBM problem used by Hydra was briefly described. The levelsymmetric LS16quadrature set was exclusively used in all calculations in this study. With respect to the spatial discretization, a Cartesian mesh grid configuration featuring 77 × 27 × 166 cells, respectively, in the radial, toroidal,and axial directions was used, yielding a total number of 345,114 mesh cells. The multigroup data library FENDL/MG-2.1 [32] and the continuous-energy data library FENDL/MC-2.1 [32] were used for Hydra and MCX calculations, respectively. The FENDL/MG library contains multigroup cross section data in the 175-group Vitamin J energy structure.The BBC and TRANSX processing codes[33] were used to prepare the 175-group P0transport-corrected cross sections for Hydra. A mono-energetic 14.1 MeV(corresponding to the eighth energy group of the 175-group structure) surface source was placed on the left surface of the DFLL-TBM calculation model. The source value was derived from the neutron wall load of 0.78 MW/m2according to Ref. [34]. All calculations of Hydra and MCX were carried out on a Linux cluster comprised of 16 computational nodes.Each node has two Intel Xeon E5620 processors and 32 GB of memory storage.

Table 1 Geometry and material configuration of the DFLL-TBM

4 Calculation results of the DFLL-TBM

The neutron flux, TBR, and shielding performance quantities including the neutron damage in terms of DPA production rates and the helium and hydrogen production rates in the steel components were evaluated using Hydra with the computational parameters presented in Sect. 3 and compared to the MCX solutions.The relative differences of the one-step method based on Hydra were also evaluated against the MCX reference solutions and compared to those obtained by the two-step approach. Note that all results were normalized to a neutron wall load of 0.78 MW/m2according to Ref. [34].

4.1 The neutronic results

The TBR is the most important quantity to measure the tritium production capability of TBM. In this study, the TBM for DFLL-TBM was calculated using the following equation:respectively. The relative difference was mere - 0.65%,demonstrating that Hydra produces the correct TBR result.A net TBR value larger than 1.0 is required to ensure tritium self-sufficiency.In fact,adequate margins in excess of unity must be attained taking into account of various tritium losses and uncertainties. A minimum TBR value of 1.2 was chosen to fulfill the tritium self-sufficiency requirement [3], taking into consideration of the nuclear cross section data uncertainties (2%), blanket ports (5%),6Li burn-up (4%), engineering design assumptions (3%),and tritium losses in the fuel cycle (5%). Clearly, the current TBM result of 0.3574 does not satisfy the tritium self-sufficiency requirement. With the assumption that ITER would be operated at a 22% duty circle (400 s burn length within a pulse of 1800 s), a total amount of 3.5 mg tritium can be produced per day in the DFLL-TBM.

Figure 8 shows the total neutron flux profile along the radial direction obtained by Hydra and MCX,respectively.Generally speaking,a good agreement is achieved between the result obtained by Hydra and the reference solution. In most regions, the relative errors were within 1%, and the largest relative error was only 3.12%, occurring in the beryllium armor. This large error was attributed to the strong anisotropy of the neutron scattering in this region which requires high-order scattering cross sections to be used, while only P0scattering cross sections were adopted in this study.

4.2 Comparison between the one-step and two-step methods

Figure 9 compares the relative differences of the total neutron flux obtained by Hydra in each material zone,with those from the two-step approach. The relative differences were evaluated against the MCX solutions. Clearly, Hydra produced results of higher accuracy compared to the twostep approach in most regions.The two-step approach used the homogenization technique to condense fine-group cross sections to coarse group ones, during which abundant information including detailed spatial and energy distribution of the neutron flux was lost. In contrast, the Hydra code directly utilized the 175-group cross sections in the calculation and hence maintained higher solution accuracy and resolution compared to the two-step approach.

4.3 Nuclear shielding results

During fusion plant operation, the high-energy fusion neutrons produced from plasma induce two damage mechanisms. The first mechanism is due to the displacement of atoms from their lattice positions as a result of collisions. The second mechanism is due to the gas production as a consequence of various nuclear reactions mainly of the (n, p) and (n, α) kind. Using the corresponding cross sections obtained from the FENDL/MG-2.1 data library and the neutron flux profile, the DPA, helium and hydrogen productions rates within structural materials were calculated assuming a 0.22 duty factor (400 s burn length within a pulse of 1800 s) within a whole year with 500 MW of thermal power. The DPA profile is illustrated in Fig. 10, and the H and He production rate results are shown in Figs. 11 and 12, respectively.

As can be seen, all neutron damage quantities reached their maximum in the FW1 region, where neutron flux spectrum is particularly hard due to the high-energy neutrons emerged from the plasma and then decreased monotonically along the radial direction by virtue of the rapid attenuation of the fast neutron flux. The maximal relative difference of these neutron damage quantities compared to the reference solutions (2%, 4.9% and 4.9%for DPA, H and He production rates, respectively) also occurred in the FW1 region. Nonetheless, in most regions the relative differences were less than 1%.In summary,all the neutronic shielding quantities obtained by Hydra agreed very well with the reference MCX solutions. This demonstrates that the discretized calculation model in Sect. 3 is sufficient to yield solutions with satisfactory accuracy.

Exposure to neutron radiation may degrade the mechanical properties of structural materials in the TBM.For example,the accumulation of helium and hydrogen can result in radiation hardening which strengthens the material while subsequently embrittling it. As a result, in order to maintain its mechanical properties, a maximal radiation limit must be set for the structural material. This is especially the case for material in the FW zone, as it directly faces the high-energy (14.1 MeV) neutrons emitted from the plasma. By using the maximum DPA limitation of 10dpa [35-37], it may be deduced that the current DFLLTBM design fulfills the operation shielding requirement of ITER for four years.

4.4 The parallel computing performance

The parallel computing performance of Hydra was also investigated using the DFLL-TBM model. In this study,two widely used performance measures were adopted,namely the speedup ratio and parallel computational efficiency (PCE). The speedup ratio is defined to be the ratio of the time to execute the computational workload W on a single processor to the time on N processors:

where τ1is the time to execute the workload on a single processor,and τNis the time to execute the same workload on N processors. The PCE is defined as:

Table 2 lists the calculational time,speedup ratio,and PCE results of Hydra for the DFLL-TBM model with different numbers of processors.

Table 2 demonstrates that computational time can be effectively reduced by employing parallel computing. The speedup ratio increased as expected with a greater number of processors being used, while the PCE results decreased due to the growing communication overhead for exchanging the boundary flux between subdomains.Nonetheless, the PCE result for the calculation with 99 processors was still 82.8%,suggesting that the Hydra code possesses good parallel scalability and can thus be applied for 3D neutronic analyses of TBMs with high computing efficiency.In comparison,it took approximately 9 h for the MCX code to accomplish the neutronic calculation for DFLL-TBM with 108particles being simulated using 16 processors.

Table 2 Parallel computing performance results of Hydra for the DFLL-TBM model

5 Parameter studies

The primary goal of the TBM program is to check and validate tritium breeding self-sufficiency capability of future breeder blanket designs for ITER. As such, understanding the tritium breeding behavior of TBM plays a vital role in promoting the development of breeder blanket designs. For this reason, parametric studies have been carried out by varying the thickness of the tritium breeding zone and the first wall beryllium layer to quantify their influence on tritium breeding performance.

5.1 Effects of the tritium breeding zone thickness

The effect of tritium breeding zone thickness on the TBR was analyzed.Seven different cases were prepared by continuously increasing breeding zone thickness from 20 to 80 cm. For each of these seven cases, the DFLL-TBM model was calculated using the Hydra code. Figure 13 shows the TBM result as a function of tritium breeding zone thickness.

Figure 13 shows that the tritium breeding performance of DFLL-TBM gradually improves as tritium breeding zone thickness continuously increases.Additionally,it also shows that the TBR converges with an increase in tritium breeding zone thickness. In fact, due to the attenuation of the neutrons, the first two zones, LiPb1 and LiPb2, contribute more to tritium production compared to the third zone, LiPb3. The TBR result will reach a maximum at a certain tritium breeding zone thickness. Beyond this point,any further increases to zone thickness will no longer improve the tritium production. In other words, the neutrons passing through the TBM have already been saturated by the tritium breeder lithium.This suggests that the use of tritium breeding zone of excessive thickness should be avoided in terms of neutron economy. In this study, a tritium breeding zone of 60 cm was selected as a starting point for the sensitivity studies of the beryllium layer thickness to tritium production.

5.2 Effects of the Be armor thickness

The (n,2n) reaction cross section of beryllium was considerably large in the high-energy groups.Accordingly,Be is usually used as the neutron breeder in many TBM designs and placed in front of the FW. In this work, the effect of the beryllium armor layer thickness on tritium production performance was investigated. Nine different calculational models were developed by continuously increasing Be armor layer thickness from 2 to 10 mm.The TBR dependence on thickness is shown in Fig. 14. As expected, thickening the Be armor layer can improve blanket TBR production performance due to the neutron breeding capability of beryllium.However,the slope of the TBR curve decreases with increased Be armor layer thickness, suggesting that a Be armor layer of excessive thickness should not be used in realistic TBM design.

6 Conclusion

In this study,the neutronic and shielding performance of the DFLL-TBM was assessed using the massively parallel 3D transport code, Hydra, with the FENDL/MG-2.1 data library. This assessment included the total flux profile along the radial direction,TBR and tritium production rate,DPA, and helium and hydrogen production rates in structural materials. The level-symmetric angular structure and the structured Cartesian mesh grid were employed to discretize the 3D DFLL-TBM problem.The calculation model was examined against the MCX code with the FENDL/MC-2.1 data library.

In most material zones, the results obtained by Hydra were in good agreement with the reference MCX solutions within 1%,indicating that the calculation model is correct,and that the Hydra code can produce accurate and faithful solutions when applied to realistic neutronic analyses for the TBM problem.A comparison of the relative differences in the results between the one-step approach based on Hydra and the two-step method suggests that the former method possesses higher accuracies as the need for a homogenization process is eliminated. The parallel computing performance of Hydra was also investigated using the DFLL-TBM model. Numerical results show that computational time can be significantly reduced through the use of parallel computing, while Hydra maintains a good parallel efficiency when using approximately 100 processors.

The neutronic calculation results show that the current DFLL-TBM design with a TBR value of 0.3574 does not meet the tritium self-sufficiency requirement.The radiation shielding results indicate that DFLL-TBM fulfills the operation shielding requirement of ITER for four years.

Lastly, parameter studies of the tritium breeding performance with respect to the Be armor and breeding zone thicknesses were conducted, providing a foundation for further improvement to the current DFLL-TBM design.