Pin-by-Pin Coupled Transient Monte Carlo Analysis Using the iMC Code

In this article, we present a coupled multi-physics Monte Carlo reactor transient analysis framework implemented in the KAIST Monte Carlo iMC code. In the multi-physics framework, the time-dependent neutron transport calculation and the transient heat transfer analysis are done based on the predictor–corrector quasi-static Monte Carlo method and the three-dimensional finite element method, respectively. Using this high-fidelity analysis framework, we demonstrated the negative temperature feedback effect in two pressurized water reactor (PWR) transient scenarios. First, a 3-D burnable absorber-loaded fuel assembly was considered with all reflective boundary conditions. In this simple problem, a positive reactivity-induced transient was analyzed to characterize the reactor responses in view of the pin-wise power and temperature distribution. Second, the iMC multi-physics analysis is applied to a control rod withdrawal transient in the TMI-1 mini core problem, and detailed time-dependent results were provided and compared with the Serpent/SUBCHANFLOW analysis. In both cases, independent MC runs were performed to quantify the uncertainty of the multi-physics MC transient analysis.


INTRODUCTION
It is a non-disputable fact that the nuclear reactor is a complex multi-physics system with neutronics, thermal hydraulics, thermo-mechanics, and chemistry interdependency. It is also well-known that a stand-alone physics simulation cannot predict the accurate and reliable dynamic behavior of a given nuclear system, yet many works had been focused on specific physics simulation mostly due to the lack of computational power.
Nowadays, a great portion of researchers in the high-fidelity reactor physics community is focusing on the establishment of a coupled multi-physics analysis framework for transient situations. The time-dependent coupled multi-physics simulations are largely based on neutronics codes with the diffusion or deterministic transport method coupled with reactor hydraulics codes (Fiorina et al., 2015;Cherezov et al., 2020;Park et al., 2020;Wang et al., 2021). In most cases, reactivity feedback effects have been deterministically applied directly to the reactivity, while the feedback constants were independently evaluated from Monte Carlo branch calculations. It is reasonable that the early multi-physics approaches are based on the deterministic neutronics methods as the time-dependent neutronics analysis itself is considered computationally exhaustive.
Despite its computation cost, the explicit Monte Carlo neutron transport method is generally accepted as the most accurate and reliable tool that takes into account various reactor feedbacks that possibly affect the dynamic behavior of a nuclear system in the most explicit way. For pursuing the high-fidelity simulation of transient nuclear systems, a few attempts have been made successfully by Serpent 2 and Tripoli-4 in combination with the SUBCHANFLOW thermal-hydraulics (TH) code (Ferraro et al., 2019;Ferraro et al., 2020a;Ferraro et al., 2020b;Ferraro et al., 2020c;Ferraro et al., 2020d). To our best knowledge, they are so far the only published research outcomes regarding the time-dependent Monte Carlo transport analysis coupled with TH features.
Unlike the Serpent/SUBCHANFLOW framework, the iMC code utilizes intra-pin power distribution for the temperature analysis to provide an enhanced fuel integrity study for a complex fuel element such as the centrally shielded burnable absorber (CSBA) (Nguyen et al., 2019). From a previous study , we found the use of detailed power distribution critical for accurate temperature analysis. We also suggested using the intrapellet power shape from a separate pin-cell transport calculation in which the effect of neighboring guide thimble can hardly affect the temperature distribution (Kim, 2022).
In this article, we present the iMC simulation results on pressurized water reactor (PWR) problems using the multiphysics coupled scheme. In Section 2, the basic framework of the coupled numerical analysis is introduced. The formulation of the predictor-corrector quasi-static Monte Carlo (PCQS-MC) for the transient neutron transport method, the threedimensional (3-D) finite element transient heat transfer analysis, and a simplified coolant model for the current analysis are explained. The numerical results are presented in Section 3, showing that the negative temperature effect of fuel and coolant on reactivity successfully moves the given reactor system to a steady state in a reactivity insertion transient.

FRAMEWORK OF TRANSIENT MULTI-PHYSICS REACTOR ANALYSIS IN IMC 2.1 Predictor-Corrector Quasi-Static Monte Carlo Method
In this section, we discuss the key features for establishing the PCQS-MC method. Guo et al. (2021) implemented the PCQS-MC framework in the RMC based on random samplings of delayed sources within a given time bin and kinetic parameter polynomial fitting method. Meanwhile, Jo and Cho (Jo et al., 2016) used analytic linear interpolation of delayed fission source and exponential transformation method. The iMC adopted the latter approach as we found it mathematically concrete and straightforward to implement.

Transient Fixed Source Iteration for Prediction Step
The time-dependent neutron transport equation includes the time derivative of flux term which is often disregarded in steady-state analyses by defining the multiplication factor. The complete governing equation in the space and time domain requires the following transport equation and precursor concentration equation: (1) Here, ψ is the neutron flux, v is the energy-dependent neutron speed, k 0 is the initial neutron multiplication factor, β is the aggregate delayed neutron fraction, C d is the precursor concentration of delayed group d, λ d is the decay constant, G d is the number of precursor group, andχ p and χ d are the prompt and delayed neutron fission energy spectrum, respectively. The operators are defined as follows: where σ t , σ s , and σ f are the macroscopic total, scattering, and fission cross section, respectively. Applying the implicit Euler scheme in the time domain gives the discretized form of transport equation: Rearranging the equation and sorting in terms of time-step flux leads to the following transport equation.
Here, the modified transport operator for PCQS formulation is defined as follows: Meanwhile, the delayed neutron source distribution at t t s can be determined by integrating Eq. 2 in the given time interval.
The integration term in Eq. 10 can be converted into an approximated form by linearly interpolating the fission source density in the time bin.
Applying Eq. 11 to Eq. 10 provides the following form of the delayed neutron source distribution: where weight factors are defined as follows: Substituting Eq. 12 to Eq. 7 results in the following modified transport equation with three delayed neutron source terms: The left-hand side of Eq. 16 represents the transport operation of the given particle, while the right-hand side terms are the sources. Eq. 16 can be rewritten to be more understandable in terms of the standard Monte Carlo source iteration framework.
where ℓ denotes the iteration step. The last two source terms with a superscript (ℓ − 1) are iteratively updated at every PCQS source iteration step. The first three source terms are sampled from the previous time step, often from the last iteration step.

Exponential Transformation
Due to the presence of delayed neutron precursors, the reactor period during transient becomes long, which makes the nuclear reactor controllable. However, this caused the simulation to cover an accordingly long period, leading to a choice between accuracy and computation time depending on the number of time steps. If one chooses a large time step size to reduce the computation time, the truncation error becomes an issue with conventional discretization schemes. The exponential transformation method applied to the reactor kinetics equation has been shown to greatly reduce truncation error caused by time discretization (Reed and Hansen, 1969;Reed and Hansen, 1970). With a properly chosen frequency in the time domain, the number of time bins for a given physical time simulation can be reasonably small since the allowable time step size is increased. Here, the change of variable is introduced for neutron flux in a time interval t ∈ [t s−1 , t s ] as follows: Then, the time derivative of the flux becomes where γ s is the exponential transformation frequency. By substituting Eq. 19 to Eq. 1 and applying the implicit Euler method in the time domain, the following discretized form of the transient fixed source transport equation at time step t s is obtained.
where Frontiers in Energy Research | www.frontiersin.org March 2022 | Volume 10 | Article 853222 Frequency γ s can be estimated by assuming that the frequencies do not change a lot over the time interval Δt. Then the approximation can be expressed as

Point Kinetics Model for Correction Step
The basic motivation of the correction step is to provide a better amplitude value calculated from a smaller time step value. The point kinetic equation can be derived by factorizing the neutron angular flux into the amplitude function n(t) and the shape function φ( r, E, Ω, t): Here, the shape function is normalized based on the initial angular flux distribution as follows: We obtain the above equation by assuming is an arbitrary weighting function.
Substituting Eq. 24 to Eqs 1, 2, and performing weighted integration over space, angle, and time results in the following point kinetic (PK) equations for the amplitude function and the weighted integrals of precursor density functions: where ρ is the reactivity, β is the delayed neutron fraction, Λ is the neutron generation time, λ d is the delayed neutron precursor decay constant, and C d is the precursor concentration. The PK parameters are defined as where W is an arbitrary weighting function, and a constant weighting function is used in the current iMC. The above PK parameters are tallied and averaged through PCQS source iterations. After the source iterations are completed, the PKE is solved based on the tallied PK parameters and PCQS micro time step δt. Here, α p (t) is linearly interpolated between time t i−1 and t i . For the discretization of the differential equation, the implicit Euler method is preferred for numerical stability.
The shape function at time step t s is calculated based on the predicted flux distribution as follows: where normalization factor Z(t s ) is defined from the normalization condition of Eq. 25.
Finally, the corrected flux distribution for the next source iteration is determined by using the amplitude function n(t s ) in the following way:

Finite Element Fuel Temperature Analysis
Based on the cell-wise power density tallied from the MC transport simulation, the FEM heat transfer calculation is performed for the fuel element temperature evaluation. The 3-D time-dependent heat conduction equation is obtained as follows: where T is the position-dependent temperature, Q is the heat source, k is the conductivity, and C p is the heat capacity.
For the finite element analysis, a linear interpolation function is applied for each tetrahedron with an interpolation function (N) defined for four nodes: where V is the tetrahedron volume and (x i , y i , z i ) are coordinates of node i. From the interpolation function, the temperature gradient in the element is obtained in terms of nodal temperatures as Using the Galerkin method, Eq. 36 is rewritten in the following form: Applying Eq. 39 to Eq. 40 gives a system of linear equations defined for every nodal temperature.
where ℓ denotes the time step index and For the case of defining one million tetrahedrons, approximately 200,000 nodal temperatures are defined. To solve such a large matrix equation, the iMC uses the Intel math kernel library dgesv routine to achieve the least burden.

Coolant Model
Before coupling with an all-inclusive subchannel program, an internal coolant model is considered for a preliminary evaluation. The active-core region coolant is lumped into a point model, considering inlet, outlet, and average temperature only. The coolant temperature and the corresponding density is the major driving factor of coolant reactivity feedback, and this is governed by the average coolant temperature in this simulation model. The lumped coolant model is illustrated in Figure 1, where the control volume envelopes the entire coolant region.
In a steady-state condition, the energy balance equation is as follows: and the outlet coolant temperature is simply obtained as where _ m is the coolant mass flow rate, c p is the specific heat capacity of the coolant, h is the coolant enthalpy, T wall is the cladding wall temperature, T f is the coolant bulk temperature, and _ Q is the thermal power of the fuel element. The average coolant temperature is defined as an average of the inlet and outlet coolant temperatures: For a time-dependent problem, the energy balance equation for the given control volume includes an additional time derivative term: Assuming a constant mass flow rate and the total mass is also constant as ρV M, the time derivative of the control volume's enthalpy is where c p ≡ zh zT f . Applying Eq. 46 to Eq. 45, the time-dependent heat balance equation becomes To numerically solve the equation, we may apply the implicit Euler scheme: The coolant temperature at the current time step is then expressed as follows:

Coupled Analysis Framework
The data exchange between the neutronics and the thermalhydraulics part in the iMC code does not require an additional file exchange protocol since each calculation module is implemented in a single platform. The neutron transport simulation generates a user-specified thermal power distribution which is to be used in the heat transfer calculation. The thermal-hydraulics module performs the heat transfer analysis to provide the temperature evolution through the time step for the next time step neutron transport simulation. Before the onset of the transient simulation, the system is needed to be in a steady state. The reactor system is first assumed to have a nominal constant temperature, and a steady-state neutron transport simulation calculates the local power distribution and global reactivity based on the arbitrary system condition. Using the neutronics output, a subsequent steady-state thermal-hydraulics simulation determines the corresponding temperature of all materials of interest. The temperature is then used in the next neutron transport simulation. Through this steady-state coupled iteration, material temperatures on a pin-by-pin basis are determined with the steady-state reactor power distribution.
Once the steady condition of the given reactor system is found, the transient simulation is commenced. Just like in the previous steady-state condition search iteration, the neutronics and the thermal-hydraulics coupled analysis is performed in every time bin, except the simulation schemes are based on the timedependent formulation. The overall calculation flow of the coupled analysis is described in Figure 2.

CSBA-Loaded Fuel Assembly
For a preliminary study on the multi-physics feedback transient analysis, a CSBA-loaded fuel assembly is designed. In this model, a 17-by-17 PWR fuel assembly is composed of CSBA-loaded fuel pellets in place of the conventional plain UO 2 fuel. Also, four control rods are initially inserted into the fuel assembly as a means to impose a reactivity transient. A discretized single CSBA pellet model and the CSBA-loaded assembly are illustrated in Figure 3.
The mesh was generated from Gmsh (Geuzaine and Remacle, 2009) software for the finite element heat transfer analysis. Detailed information about the CSBA model is presented in the study by Kim and Kim (2021b).
To simplify the problem, a single layer of CSBA-inserted fuel is considered and the assembly is considered to be subjected to all reflective boundary conditions in four radial sides, top, and bottom. In the iMC, a user can provide reactivity transient scenarios in two different ways: changing material cross section by mixing or replacing with predefined materials or deforming the geometry by tuning surface parameters. In either case, the time-dependent parameter values are given in a user-defined function of time. To simulate the control rod withdrawal transient in this problem, we linearly mixed the coolant with  the control rod material instead of geometric movement. The time-dependent mixing ratio is given in terms of volume fraction, and the corresponding new density is calculated in the iMC. Before the transient simulation, a series of feedback iterations is performed to find an initial steady-state reactor system with the initial thermal power of 36 W/gU. The initial intra-pin temperature distribution and coolant average temperature are set as a constant value for the first steadystate transport simulation. With the obtained power distribution, the detailed temperature is updated. The iteration quickly converged to a steady-state after 2-3 iterations. The obtained pin-power distribution and pin average fuel temperature distribution are shown in Figure 4. The CSBA model used in this simulation is described in Table 1, and material properties used in the analysis are identical to those in the study by Kim and Kim (2020).
Starting from the initial material temperature and coolant density, the PCQS-assisted Monte Carlo transient calculation is done. The control rod material is linearly mixed with water from 0.5 to 1.5 s up to 20% of water volume fraction; 25 independent batch runs were performed to evaluate the uncertainty of the simulation. The reactor thermal power profile over the transient is shown in Figure 5. Without adequate temperature feedback from fuel, burnable absorber, and coolant, the reactor power increases exponentially with the given positive reactivity as the dashed line in Figure 5. However, the negative feedback effects of fuel and moderator temperatures suppressed the power excursion right after the end of additional reactivity insertion; the reactor power started to converge toward another stable state. Figure 6 presents the fuel and moderator average temperature evolution with respect to the given thermal power transient. The fuel temperature response is rather prompt, although there is a  slight lag to reach a stable state. This lag is largely due to the thermal inertia of the fuel and the heat dissipation delay from the fuel to the coolant. The coolant temperature response was clearly slower than that of the fuel since it is the last material of the heat transfer process. These temperature responses are also shown in local power and temperature in Figures 7, 8. The peak pin-power occurred at 1.5 s and decreased afterward, while the pin-wise fuel temperature monotonically increases to reach a plateau. Table 2 shows the computation burden for major functions in a fraction of total computation time. The particle transport occupies most of the resources in the multi-physics simulation. The source list setting and combing takes less than 1% of the total burden, and the PK equation solving time is negligible. The FEM heat transfer calculation for each fuel pin is rather significant, implying a possibility of acceleration and algorithmic improvement for the matrix operation for larger problems.

PWR Mini-Core Problem
This numerical study is for the verification of the iMC multiphysics coupled transient analysis framework. The only currently available Monte Carlo coupled analysis result is provided by the Serpent/SUBCHANFLOW framework. The TMI-1 mini-core problem provided in Ref (Ivanov et al., 2013) is considered for the iMC analysis, and the results are compared against the Serpent/SUBCHANFLOW analysis result presented in Refs. (Ferraro et al., 2019) and (Ferraro et al., 2020a). The TMI-1 mini-core problem consists of 15-by-15 PWR fuel assemblies in a 3-by-3 configuration as shown in Figure 9. In the center fuel assembly, 16 Ag-In-Cd control rods are initially fully inserted.
The coolant contains soluble boron to make the system critical. The critical boron concentration (CBC) depends on the thermal-hydraulic model since the temperature and grid structure of the coolant channel are directly related to the borated coolant density. In the study by Ferraro et al. (2019), Serpent 2 and Tripoli-4 used a point-wise coolant model and neglected time delay in heat transfer from the fuel to the coolant, and the iMC obtained a similar value as in Table 3. When the coolant model is based on the detailed subchannel analysis code, the CBC value is lowered by about 200 ppm from the simple model due to the aforementioned reason.
The simplified TH model in Ref. (Ferraro et al., 2019) assumes no thermal inertia of the fuel element, exhibiting no time delay of heat deposition inside the fuel by neglecting the time derivative term of the transient heat transfer equation. Since the iMC considers time-dependent heat transfer inside the fuel element without any approximation (but point model for the coolant), comparing the iMC analysis result with the simplified fuel heat transfer model is not adequate. Instead, we proposed to compare with the Serpent/SUBCHANFLOW coupled analysis model while the amount of static reactivity insertion is set to be identical. The static reactivity inserted in the Serpent/SUBCHANFLOW analysis is 354 pcm, with control rods withdrawn by 30 cm. We evaluated the equivalent control rod withdrawal length in the iMC that provides a similar static reactivity ( Table 4). The critical boron concentration of the iMC model is determined on a trial basis iteration.
The reactor transient starts from 0.2 s by pulling out the 16 control rods from the active core with a constant speed precalculated in Table 4 for 1 s. The reactor is axially divided into  Figure 10 and Table 5, respectively. As shown in Figure 10, without the negative temperature feedback of the reactor system, the thermal power exhibited an exponential excursion with the positive reactivity insertion. The iMC and Serpent simulation agreed well with each other except for the power evolution in response to the initial control rod withdrawal. This is because the control rod differential worth is slightly different in the two cases, albeit identical total reactivity insertion. Since the differential control rod reactivity worth becomes more significant as it approaches the active core center, the power growth is steeper in the iMC model as the withdrawal length is longer. The total reactivity insertion was similar in both cases (not exactly the same due to temperature difference), so the peak thermal power is FIGURE 10 | Comparison of Serpent/SUBCHANFLOW and iMC multiphysics using TMI-1 mini-core problem.  evaluated to be very similar in the two codes. Figure 11 illustrates the local power transient of the mini-core problem. After 0.2 s of stable null transient, the control rod assembly is withdrawn from the bottom of the core for 1 s at a constant speed and stopped afterward. The axial power peak shifts toward the bottom of the core, and the local power transient behavior well agrees with the point reactor power result. Meanwhile, the time-dependent fuel pin temperature shown in Figure 12 monotonically increases during the period of the entire transient similar to the CSBA-loaded fuel assembly study. Also, the temperature peak skews toward the bottom of the core as the upward withdrawal of the control rod bank induces more bottom-skewed axial power profile.
The computation time required for the simulation is relatively large compared to the steady-state simulation even for the Serpent. The Serpent Monte Carlo transport simulation is based on the dynamic Monte Carlo method, so the computation time is less affected by the time step size. Although the two analyses were performed in different calculation conditions, the iMC transport simulation time is still much longer than that of the Serpent using a longer macro time step size. However, this performance gap can be reduced by adopting the coarse-mesh acceleration methods to reduce inactive cycles and the ray-tracing improvement of the iMC transport function itself. Future effort on algorithmic enhancement is needed for such a demanding process along with a sufficient amount of computation

CONCLUSION
In this work, we presented a multi-physics framework implemented in the iMC. The PCQS-MC-based transient neutron transport simulation and the finite element heat transfer are solved in a coupled framework to account for the temperature feedback effects in the time domain. The multi-physics framework is tested for a CSBA-loaded PWR fuel assembly when external reactivity is inserted via control rod removal from the system. The feedback correction of material cross-sections with adjusted temperature and density suppressed the additional reactivity and led the system to a stable state. The temperature responses of the system showed a slight lag from the initial perturbation due to the heat transfer delay, resulting in the power overshoot before the stabilization. The iMC multi-physics framework was also compared with the Serpent/SUBCHANFLOW coupled analysis framework on the TMI-1 mini-core problem for verification. The calculated time-dependent power evolutions from the two codes matched well, except for a minor discrepancy in the initial reactivity insertion period due to the different differential worth of the control rod. The iMC analysis framework is especially useful for the 3-D complex fuel element as it can utilize the Gmsh-generated unstructured mesh grid for spatial discretization. For a more realistic reactor system simulation, the iMC multi-physics framework will adopt a subchannel analysis program in the near future. As we have demonstrated, such detailed analysis is essential for the accurate safety analysis of various fuel designs. We expect the applicability of the iMC multi-physics feature can be far extended for advanced PWR fuel elements that are continuously being developed and even for the molten salt reactor analysis.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HK developed the program and carried out the result under the supervision of YK.