Neutronics analysis code MagicMC

MagicMC is a neutron/photon transport calculation code based on the MC method, which can complete the modeling of a variety of complex geometries. It has developed a special time-dependent particle transport module for transient neutronics/thermal-hydraulics coupling calculations.

MagicMC coupled with the on-the-fly cross-section generation method can quickly and accurately evaluate the impact of the Doppler effect [25] on the resonance cross sections. Owing to the TRIGA reactor’s pulse-transient operation status, the reactor changes from critical to supercritical and then to subcritical and finally reverts to its critical condition within an extremely short duration. To realize an accurate description of the pulse process induced by the continuous insertion or withdrawal of the pulse rod, MagicMC not only developed a time-dependent geometric movement function (moving geometry method) to describe the rod movement in the reactor, but also calculated the time-dependent power changes under supercritical status through efficient control of the particle population with the combined method. The functional modules of the moving geometry, time-dependent particle transport, on-the-fly cross-section generation, and time-dependent tally developed by MagicMC for neutronic/thermal-hydraulic coupling are shown in Fig. 2.

Fig. 2Full-screen View

(Color online) The functions modules developed by MagicMC in response to neutronics/thermal-hydraulics transient coupling

2.1.1

Moving geometry method

In this study, to realize MC calculations with moving geometries usually occurring under TRIGA reactor transient conditions, MagicMC decomposes the MC model into fixed and moving geometry parts. When the particle enters the moving geometry region, the positions of the moving geometries are recalculated based on the current particle time variable and moving geometry motion parameters (i.e., translation and rotation), and the particle transport calculation is restarted. At the end of the transport time step, the particle time variable and positions of the moving geometries were recalculated, the cell parameters were updated, and the next transport step was started. This cycle is repeated until the particles disappear or leave the moving geometry region. Figure 3 illustrates the particle transport process in a moving geometry.

Fig. 3Full-screen View

(Color online) MC particles transport process with moving geometry method

2.1.2

Time-dependent particle transport with the combing method

In this study, the combination method is introduced as a population control algorithm, which realizes the control of the number of particles by combining the particle weights when the neutron population is increased drastically in the TRIGA pulse-transient operation status. Assuming a total of K particles in the system at moment t, the total weight of the particles accumulated by the system is given by(1)A data set was generated according to the following method:(2)where ξ is a uniformly distributed random number from 0 to 1 and N is the number of particles to be maintained.

When cm falls into the weight interval wi, the corresponding particle is copied and assigned a new weight value:(3)Figure 4 shows the process of compressing K particles into N particles with weight . Meanwhile, by defining the integer j, particles with weights wi may produce j particles or j+1 particles after combining(4)The probability of producing j and j+1 particles is(5)(6)Then, the expected values of the particle weights before and after the compression of the particle with weight are(7)Therefore, the weights of the particles compressed using the combined method are conserved in the statistics.

Fig. 4Full-screen View

(Color online) The combing method is used to compressing K particles into N particles

Subchannel analysis code SubChanFlow

SubChanFlow(SCF) is a subchannel-level thermal code for the steady-state and transient analyses of a reactor developed at KIT [26]. This is written entirely in Fortran [27]. Mass conservation, energy conservation, axial momentum conservation, and transverse momentum conservation equations are used to construct the solution model. At present, SCF can define a variety of coolant fluids such as water, lead, lead-bismuth, sodium, etc. Simultaneously, based on the subchannel method, SCF can flexibly divide the control body in the reactor. In theory, SCF can complete the model construction of all rod-fuel reactors. Figure 5 shows the input information and computational flow of SCF.

Fig. 5Full-screen View

(Color online) The input information and the computational flow of SCF

Neutronics/thermal-hydraulics coupling method

2.3.1

Steady-state coupling method

Because of the absence of time-dependent treatment methods for the steady-state problem, the steady-state coupling method based on MagicMC and SCF differs significantly from the transient coupling method. The steady-state coupling process uses an overrelaxed iterative method. Steady-state coupling starts with an initial thermal-hydraulic calculation with the power set to a uniform distribution. The calculated temperature and coolant density fields were then passed to the input model of the MagicMC code, and the fuel temperature and coolant density were reset according to the mapping relationship, as shown in Eq. 8. The MagicMC code performs neutron transport to calculate the power distribution and passes it to the SCF model to correct it. Thus, the iterative calculation outputs the results when the convergence conditions are satisfied. The steady-state MagicMC-SCF coupling method is illustrated in Fig. 6. This method uses the infinite norm of the fuel temperature, calculated twice using adjacent iterations, as the convergence condition [28, 29], as shown in Eq. 9.(8)(9)where represents the temperature of the reactor fuel, ω is the relaxation factor, ranging from 1 to 2, is the calculated fuel temperature for SCF, n indicates the number of iterations, i,j,k represent the position information of the fuel element, is the relative error calculated at the nth iteration, ε is the convergence criterion. In the MC calculations, the number of simulated particles strongly influenced the convergence of the results. If the number of particles in the calculation is not controlled, it can lead to the method falling into a situation where it cannot converge when the convergence conditions are set too harshly. In addition, considering the error in the initial assumption of uniformly distributed power, an adaptive population control strategy was developed to adjust the number of neutrons in the MC calculations to satisfy the convergence condition in neutronics.(10)where Nn is the number of particles simulated in the n th iteration of the computational simulation, n_max is the maximum number of iterations used to eliminate errors caused by uniformly distributed power.

Fig. 6Full-screen View

(Color online) Steady-state coupling method based on MagicMC-SCF

2.3.2

Transient coupling method

In the transient analysis, because the time-dependent changes of the neutronics and thermal-hydraulics make the coupling process relatively complex, semi-implicit coupling is used as the time-coupling strategy for MagicMC-SCF coupling in this work. The transient-coupling method is illustrated in Fig. 7.(11)(12)where t is a time point; Δt is the time step; is the infinite norm of the temperature field between adjacent time steps; is the direct result of the coupling method in the n’th iteration of the calculation at the moment t; a is the time step increase factor; b is the time step shortening factor; N_max is the maximum number of iterations, ΔT_min and ΔT_max are semi-implicit iterative judgment conditions. As time advances, the code first calculates the infinite norm of the temperature field between adjacent time steps, and determines the convergence of the data in a single time step through semi-implicit coupling. The specific method is as follows.

Fig. 7Full-screen View

(Color online) Transient coupling method based on MagicMC-SCF

a)

When , the data between adjacent time steps are considered to have good convergence and can be output directly. This method carries out transient coupling to advance and lengthen the time step to reduce the number of data exchanges.

b)

When , the coupling method directly considers that the coupling result has not converged, reduces the time step, and returns to the calculation.

c)

When , the coupling data is considered to be between "convergence" and "non-convergence,” and the over-relaxation iteration will be carried out in the time step at this time. During the iteration of this time step, the method still performs computation and judgment of . If the conditions in (a) or (b) are satisfied, the method prioritizes the execution of the time-step computation strategy in a) or b). If the conditions in (a) and (b) are not satisfied, iterations are performed based on the temperatures computed in Eq. 11 until convergence, after which time advancement is performed.

Because the transient coupling calculation is a process advancing with time, in practice, the neutron source calculated at different time steps will theoretically change. The coupling method updates the neutron source in each time step to minimize its influence on the coupling results.

Neutronics model

Two different core arrangements are selected as numerical model for steady and transient states. Core 134 is selected as the steady-state model, and core 231 is selected as the research object for the transient model. The core structure diagram constructed using MagicMC is shown in Fig. 8(e). In core 134, all control rods except the regulating rod are fully withdrawn. In core 231, both the safety and shim rods are fully withdrawn, the transient rod (i.e., pulse rod) is fully inserted into the core, and the regulating rod is partially withdrawn to make the core exactly critical.

For the neutronics model, to smooth the temperature and power of the neutronics model, a zirconium rod is designed at the center of the fuel element, and small sections of graphite are designed at the upper and lower ends of the fuel element. The reactor core is reactively controlled by four control rods: a safety rod, shim rod, regulation rod, and transient rod. Except for the transient rod, all the other control rods are designed with a fueled-follower. Figure 8(d) shows the design scheme of the fueled-follower control rod and fuel element. Cores 134 and 231 are designed with 12 wt% uranium and U235 with 20% enrichment fuel elements. Table 1 lists the materials and dimensions used in the neutronics model.

Thermal-hydraulics model

The subchannel division method generally uses the coolant channel as the center [33]. In 2021, Manuel Garcia of KIT analyzed the two partitioning methods centered on the coolant channel and fuel element using the VVER benchmark and concluded that negligible difference existed between the results of the two partitioning methods [29]. Because the TRIGA Mark II reactor uses ring-filled fuel elements for the reactor design, this study adopts a subchannel division method centered on the fuel elements. The core sub-channel division is illustrated in Fig. 8(f).

The TRIGA reactor uses natural circulation is used to cool the TRIGA reactor. The thermodynamic parameters of fuel U-ZrH are shown in Eq. 13 and Eq. 14 and the thermodynamic parameters of stainless steel are given by Eq. 15 and Eq. 16 [18, 34].(13)(14)(15)(16)where Cv is the heat capacity per unit volume, J/(m³·K)is the unit, and k is the thermal conductivity in units of W/(m·K).

Spatial mapping method

Spatial mapping methods can be divided into volume weighting and one-to-one mapping methods. Because no power is released from the coolant during the actual calculation, the data processing for the fuel element and coolant differ. The main data used for the mapping are power, temperature, and density. The following method is used to map the data:

Fig. 9Full-screen View

(Color online) Spatial mapping of temperature data(10 nodes in axial direction)

Steady-state coupling results

The steady-state calculation used core 134 as the research analysis model, and two different power levels, 125 and 250 kW, were calculated. The two power levels were regulated by the position of the regulating rod, and the remaining control rods were fully withdrawn. At 125 and 250 kW, the regulating rod was inserted up to 19.4 and 13.7 cm, respectively; the coolant inlet temperature was 21.6 ℃ at both the power levels. Table 2 lists the information on k_eff for the steady-state coupling of the two power levels.

Figure 10(a-b) are the pin-wise distribution of radial power factor. The standard deviation of all power factors did not exceed 0.003. It can be observed that the power factor of ring B is higher than that of ring A for both power levels, and this result follows the value of the fuel elements of rings A and B in [31]. Because the steady-state power is controlled using only the insertion behavior of the rod R, the power factors of the rod R and its nearby elements differ significantly at the two power levels. In addition, the pinwise power peak factors at 125 kW and 250 kW are 1.3734 and 1.3403, respectively, which verifies that the deeper the R-rod is inserted, the more uneven the core power distribution. To demonstrate high-fidelity power distribution, Fig. 11(a-b) show a finer radial power density distribution in the active zone of the core at 125 kW and 250 kW using an x-y grid of 500×500.

Fig. 10Full-screen View

(Color online) Radial pin-wise power factor distribution(R is regulating rod, C is shim rod, S is safety rod, and T is transient rod).(a) Pin-wise power factor distribution at 125 kW. (b) Pin-wise power factor distribution at 250 kW

Fig. 11Full-screen View

(Color online) Steady-state coupling results. (a) Radial power density distribution at 125 kW. (b) Radial power density distribution at 250 kW. (c) Normalized axial power distribution (d) Fuel temperature in different rings at 125 kW and 250 kW

The normalized power distribution is calculated in the axial direction of the reactor and the results are shown in Fig. 11(c)). The relative axial position of only the active zone of the core is calculated, with the bottom corresponding to the relative position 0 and the top corresponding to the relative position 1. The axial power peaking factors at 125 and 250 kW are 1.288687 and 1.300017, respectively.

Figure 11(d) shows a comparison of the calculated and experimental values of the fuel temperature for different rings. Table 3 compares the calculated and experimental temperatures at certain fuel locations. These temperatures are located ~7 mm from the center of the fuel element in the radial direction and their axial position is at half the height of the fuel meat. The accuracy of absolute temperature measurements in the experiments, due to calibrating uncertainty, is estimated ±3℃ [31]. The temperatures calculated in this study are consistent with the experimental values. The main error arises from the necessary model simplification required for the code-based computation.

Pulse transient coupling results

The transient calculations use core 231 as the analysis model. During the transient calculation, the safety and shim rods are fully withdrawn, and reactivity is introduced by withdrawing the transient rod for a short time (~ 80 ms). The integral value of the transient rod under core 231 is calculated to select the correct position for the transient rod. The reactivity worth of the transient rod is calculated using the positive-period method [35], as shown in Eq. 17 and the error propagated during the calculation is given by Eq. 18 [36].(17)(18)where denotes the effective multiplication factor. Subsequently, one of the control rods is withdrawn at a certain position by calculating a new k_eff. and are the fractional statistical error estimates for and k_eff, respectively.

The integral worth of the transient rod in core 231 is shown in Fig. 12(a). Compared to the reactivity value of the transient rod measured in the experiment [6], the integral value of the transient rod is roughly consistent with the experimental value, and the main error arises from the uncertainty of the model simplification construction and nuclear data.

Fig. 12Full-screen View

(Color online) Pulse transient coupling results (a) Comparison of MagicMC and experimental transient rod integral worth in core 231. (b) Neutron energy group distribution at the central guide tube at peak power for two transient conditions. (c) Variation of power over time in 2.75$ reactivity insertion (d) Variation of power over time in 3$ reactivity insertion

Based on the calculated integral worth of the core 231 transient rod, this study used a coupled code to perform a reactivity insertion pulse transient simulation. The reactor started with an initial power of 100 W and coolant temperature of 21.6 ℃. The length required for the transient rod to be withdrawn was calculated based on the integral value of the transient rod, as shown in Fig. 12(a). The fixed-source mode was used to calculate the power change in the reactor in a single time step. The shortest timescale of the power change in a single step was 0.002 ms (the minimum scale of the solving step).

The pulse power as a function of time is shown in Fig. 12(c-d). To analyze the error source, the peak power of the power curve, calculated using the coupling code, was aligned with the experimental value. Table 4 compares the calculated and experimental values of the key parameters of the pulse-coupling model. The experimental data used in this study were obtained from the JSI public TRIGA pulse experiment, which can be found on the webpage http://trigapulse.ijs.si/. By comparing the power curve and pulse parameters, the coupling method proposed in this work was found to accurately simulate the pulse transient of the TRIGA reactor. The calculated and experimental power change curves are generally consistent. Small differences exist only in regions where the power varies relatively slowly because the presence of delayed neutrons in these regions affects the power variations. In the two simulations, the relative errors between the calculated and experimental values of the peak power were 0.705 and 1.834%, respectively, and the relative error of the maximum fuel temperature did not exceed 4.5%.

Figure 13(a-b) display the temperature distribution of at different time steps under two transient coupling calculations. The temperature decreased from the inside to the outside. However, the positions of the center and transient control rod do not have a fission power release; therefore, heat cannot be transferred to these locations in a short pulse time. Consequently, their temperatures were significantly lower than those of other parts of the core.

Fig. 13Full-screen View

(Color online) (a) Element temperature distribution at different time step in 2.75$ reactivity insertion. (b)Element temperature distribution at different time step in 3$ reactivity insertion

Because the coupling code uses MC for 3D neutronics calculations, it can output neutron flux distributions for different energy groups at different time steps. The central guide tube was selected as the analysis model, and the WIMS 69-group and WIMS 172-group structure [37] were used as the energy group division method. The flux distribution per lethargy of groups 69 and 172 at the central guide tube when the peak power was reached under the two transient conditions is shown in Fig. 12(b)). The fluxes per lethargy for the different energy-group structures were in good agreement. Because the two-pulse calculations had little effect on the core structure, the relative flux distributions under different energy groups were basically the same. A distinct difference is observed only in the flux magnitude per lethargy mainly because of the difference in the peak power owing to the magnitude of the reactivity insertion.

Computing environment

All calculations were performed on a computer with an Intel Core i5-13600KF CPU and 48G GB RAM. The calculations were performed using Linux with the CPU operating at a stable frequency of 4.6 GHz. An OpenMP parallel programming model with 20 threads was used for the computations.

In the steady-state coupling, the calculations were performed with 100 inactive cycles followed by 1000 active cycles, with 5000 neutrons per generation. The convergence condition () for the coupled system is set to 0.002. For 125kW operation, the coupled system converges after 59 iterations, and the number of neutrons computed per generation increases to 22533. For 250kW operation, the convergence condition is reached after 49 iterations, and the number of neutrons computed per generation increases to 26629. The runtime for both the calculations was approximately 3 h.

In transient coupling, MagicMC used 1.5 million particles in each time bin, and the time bins under a single coupling time step were evenly divided into 100. In the adaptive time-step strategy, we set ΔT_max and ΔT_min to 6° and 2° respectively, which controlled the time-step behavior of the coupled system. Semi-implicit calculations at a single time step were considered to converge when was less than two ℃. The calculated runtime was approximately 28 h with 2.75$ reactive insertion and approximately 30h with 3$ reactive insertion, and the minimum time step for coupling was 0.0004s.