A Multi-Physics Adaptive Time Step Coupling Algorithm for Light-Water Reactor Core Transient and Accident Simulation

Abstract

A new reactor core multi-physics system addresses the pellet-to-cladding heat transfer modeling to improve full-core operational transient and accident simulation used for assessment of reactor core nuclear safety. The rigorous modeling of the heat transfer phenomena involves strong interaction between neutron kinetics, thermal-hydraulics and nuclear fuel performance, as well as consideration of the pellet-to-cladding mechanical contact leading to dramatic increase in the gap thermal conductance coefficient. In contrast to core depletion where parameters smoothly depend on fuel burn-up, the core transient is driven by stiff equation associated with rapid variation in the solution and vulnerable to numerical instability for large time step sizes. Therefore, the coupling algorithm dedicated for multi-physics transient must implement adaptive time step and restart capability to achieve prescribed tolerance and to maintain stability of numerical simulation. This requirement is met in the MPCORE (Multi-Physics Core) multi-physics system employing external loose coupling approach to facilitate the coupling procedure due to little modification of constituent modules and due to high transparency of coupling interfaces. The paper investigates the coupling algorithm performance and evaluates the pellet-to-cladding heat transfer effect for the rod ejection accident of a light water reactor core benchmark.

1. Introduction

Multiphysics analysis attracts a lot of attention of nuclear engineers worldwide as it helps to produce more realistic results in terms of reactor core safety margins. Nowadays there are many codes progressing towards the high-fidelity multiphysics simulation of light water reactors [1,2]. Recent works [3,4] demonstrate the coupling experience of MCS [5], CTF [6] and FRAPCON [7] codes for reactor core pin-by-pin analysis taking into consideration the pellet-to-cladding transfer and coolant cross-flow phenomena. To continue this work, the Computational Reactor Physics and Experiment Laboratory of Ulsan National Institute of Science and Technology (UNIST CORE lab) has developed MPCORE (Multi-Physics CORE) code [8] capable for both depletion and transient analysis of a light-water reactor core with the dynamic pellet-to-cladding heat transfer coefficient taking into account fuel rod irradiation, thermal-mechanics and corrosion phenomena. The new multiphysics system implements the external loose coupling of pin-by-pin neutronics, thermal-hydraulics and fuel performance using the external adaptive time step selection algorithm to achieve the reactor core high-fidelity simulation with prescribed tolerance. The processing of the multi-physics data producing by the MPCORE code has been discussed previously [9]. This work is devoted to the application of the adaptive time step integration algorithm, analysis of the code computational performance and evaluation of the pellet-to-cladding heat transfer effect on the safety parameters using for reactor core rod ejection accident (REA) transient [10,11,12].

MPCORE performs numerical integration of reactor core codes using Implicit Euler step-doubling method [13,14]. The time step selection algorithm estimates time discretization error at the current time step and then chooses the next step size in such a way that a given tolerance is not exceeded. For a chosen step size one transient step is calculated with this step size and compared to the result of two transient steps with half the step size. This approach is known as the two-level step-doubling method described in the works [13]. We chose this method as being more robust compared to embedded pair methods such as considered in [14,15]. The implicit equation within a time step is solved by the damped Picard method [14] with Gauss–Seidel acceleration for better convergence.

The current MPCORE version comprises three computational modules: three-dimensional two-group nodal diffusion code RAST-K [16] augmented a pin-by-pin reconstruction module, homogeneous one-phase Coolant Thermal-Hydraulics One Dimensional code CTH1D and Fuel Rod Analysis Programming Interface FRAPI designed for the external loose coupling of fuel rod analysis codes FRAPCON [7] and FRAPTRAN [17]. The description and verification of FRAPI code was presented recently [18].

The paper is organized in five sections. Section 1 is the introduction part. Section 2 describes the pellet-to-cladding heat transfer phenomena, multi-physics system modules and the multi-physics coupling algorithm. Section 3 is devoted to the analysis of the rod ejection accident (REA) from hot zero power (HZP) core, effect of the gap heat transfer coefficient (HTC) on the reactor core safety parameters, and simulation convergence and performance. Section 4 contains the conclusion and future work plan.

2. Methodology

2.1. Pellet-to-Cladding Heat Transfer Phenomena

The temperature distribution in a fuel pellet, cladding and coolant is of primary importance of multi-physics simulation because it has a great influence on the distribution of neutrons due to the Doppler effect and coolant density change. The heat transfer phenomena inside a fuel rod comprises heat conduction in the fuel pellets and cladding, and the heat transfer across the pellet-cladding gap filled with a gas under high pressure [19]. The pellet heat conduction mainly depends on the heat capacity and thermal conductivity coefficients which are the functions of local fuel burnup, temperature and density. Notwithstanding the influence of the restructuring effects such as fuel cracking, grain growth and high burnup structure formation, the effective fuel Doppler temperature used for neutronics feedback is predicted with enough high accuracy for engineering calculation using the heuristically obtained correlations. In contrary, the simulation of the heat transfer across a gap is a more challenging problem as the pellet-cladding HTC primarily depends on the radial distance between a pellet and cladding surfaces, and, consequently, affected by the thermal-mechanics and fission gas release. Therefore, the rigorous simulation of pellet-to-cladding heat transfer involves synthesis of neutronics, coolant thermal-hydraulics and fuel performance codes.

The pellet-to-cladding heat flux $q_F$ in watts per unit area is calculated as follows

$$q_F=h_{gap}(T_{fs}-T_{ci}),$$

(1)

where $T_{fs}$ is the outer surface fuel pellet temperature in K, $T_{ci}$ is the inner surface cladding temperature in K and $h_{gap}$ in W/m${}^2$K is the heat transfer coefficient between fuel and cladding. Overall, the fuel-cladding gap conductance model [20] consists of four parallel conduction routes:

$$h_{gap}=h_{solid}+h_{gas}+h_{rad}+h_{con},$$

(2)

where $h_{gas}$ is the gas heat conductance, $h_{rad}$ is the conductance by radiation from fuel outer surface to cladding inner surface, $h_{solid}$ is the conductance by solid pellet-cladding contact and $h_{con}$ is the conductance by gas convection which can be neglected as the gas is quiescent and fuel rod radial gap size is very small, about 60 $\mathsf{\mu }$m at the hot state.

The radiation heat conductance term $h_{rad}$ is described by the Stefan–Boltzmann law

$$h_{rad}=\sigma ϵ(T_{fs}^4-T_{ci}^4)/(T_{fs}-T_{ci}),$$

(3)

where $\sigma$ is the Stefan–Boltzmann constant, $ϵ$ is the emissivity factor. The thermal conduction by radiation is normally less than 1% during normal operation condition due to the limited surface temperatures. However, if cladding ballooning occurs, the distance between pellet and cladding surfaces becomes too large and radiation heat transfer can outperform the heat transfer by gas conductance.

The pellet-cladding gap filled with the as-fabricated helium gas acquiring a large thermal conductivity coefficient and gaseous fission products, primarily krypton and xenon, released from uranium dioxide during irradiation. The gas heat conductance is given by the model of Ross and Stoute [21],

$$h_{gas}=\frac{{\lambda }_{gas}}{{\delta }_{gap}+{\delta }_f+{\delta }_c+{\delta }_{jump}}$$

(4)

where the thermal conductivity of a multicomponent gas ${\lambda }_{gas}$ is divided by sum of the gap size ${\delta }_{gap}$, fuel and cladding mean roughness ${\delta }_f$ and ${\delta }_c$, and combined fuel and cladding gas temperature jump distance ${\delta }_{jump}$.

The temperature jump distance, also known as the gas extrapolation length, accounts for the temperature discontinuity near the fuel and cladding walls caused by the accommodation effect. The accommodation effect appears to be important if the gap width is small compared with the gas mean free path and gas molecules travel back and forth between walls rarely colliding with each other. For example, at 20 MPa gas pressure and 900 K gas temperature the jump distance for helium is about 10 $\mathsf{\mu }$m, for xenon about 1 $\mathsf{\mu }$m. According to the Maxwell gas–surface interaction model, gas molecules undergo the diffusive reflection associated with the complete thermal accommodation to a wall surface temperature and the specular reflection that does not change the particle’s energy. The specular reflection and the inability of gas molecules leaving a surface to fully exchange their energy with neighboring gas molecules results in a steep temperature gradient near the fuel and cladding surfaces which is computationally adjusted by ${\delta }_{jump}$ in Equation (4). The detailed gas–surface interaction model, including the fraction of specular reflection events, for different gas and wall materials is the subject of many experimental measurements and theoretical studies [22].

Thermal expansion, cladding creep under external pressure, irradiation effects, fuel relocation due to cracking all contribute to the mechanical deformation of a fuel pellet and cladding, which results in a gradual decrease in the gap width ${\delta }_{gap}$ over fuel burn-up until it disappears after 15–20 GWD/MTU [23]. The pellet-cladding mechanical interaction at lower burnups is mainly observed for reactivity initiated accident events accompanied with high temperatures owing to the differential thermal expansion between fuel pellet and cladding. When the fuel pellets come into contact with the cladding, the heat transfer coefficient $h_{gap}$ dramatically increases owing to thermal conduction through contact spots separated by microscopic valleys of the contact area. The solid conductance term $h_{solid}$ depends on the thermal conductivity of cladding and pellet, interfacial pressure and surface roughness [17,21]. Thus, owing to the presence of surface asperities, the heat transferred across the contact interface consists of three components $h_{solid}$, $h_{gas}$ and $h_{rad}$ with the primary contribution arising from the first term.

2.2. MPCORE Modules

MPCORE is designed in such a way that it is able to link any number of external codes using the interface functions demanded from each external module. The current version supports three computational modules:

• RAST-K [16] is a three-dimensional two-group nodal diffusion code used for depletion and transient analysis of light water reactor cores having square assemblies. The nodal solver is accompanied with pin-by-pin power reconstruction, micro-depletion analysis and xenon/samarium build-up modules. The two-group cross-section library and pin-wise power shape functions for the reconstruction module are computed by lattice code STREAM. This allows the fuel assemblies with asymmetric burnable rod position to be assessed [24].
• FRAPI [18] is the fuel rod analysis program interface designed for the loose coupling of the fuel rod analysis codes FRAPCON and FRAPTRAN [17] with other reactor core codes. The interface supports most of the original code options beside the two-dimensional finite element analysis of fuel rod cladding. The diagram of the phenomena simulated by the codes is presented in Figure 1.
• CTH1D is the homogeneous one-phase code for transient and steady-state simulation of coolant flow in a fuel rod channel. The code implements one-dimensional energy and mass conservation equations with the steam-water properties provided by FREESTEAM which is an open source implementation of international-standard IAPWS-IF97. In addition, CTH1D is able to simulate coolant boiling effect based on the correlations of critical heat flux and heat transfer coefficient adopted from FRAPTRAN.

The feedback exchange between computational modules is implemented according to the scheme presented in

Table 1. Exchange of feedback parameters.

From/To	RAST-K	FRAPI	CTH1D
RAST-K	-	$q_l$	-
FRAPI	$T_{fuel}$, $T_{clad}$	-	$T_{film}$
CTH1D	$T_{cool}$, $D_{cool}$	$H_{film}$	-

. In detail, RAST-K calculates pin-by-pin linear power distribution $q_l$ which is forwarding to FRAPI. The FRAPI code calculates fuel pellet Doppler temperature $T_{fuel}$, average cladding temperature $T_{clad}$ and cladding oxide film temperature $T_{film}$, the last two parameters are used as the boundary conditions for coolant thermal-hydraulics. CTH1D calculates bulk coolant temperature $T_{cool}$, bulk coolant density $D_{cool}$ and film heat transfer coefficient $H_{film}$ adopting by fuel rod analysis code FRAPI. Updated $T_{fuel}$, $T_{cool}$, $D_{cool}$ and $T_{clad}$ variables are used by RAST-K for the two-group cross-section calculation. The iteration repeats until the convergence of all the feedback parameters.

2.3. Integration Scheme

MPCORE code performs the external loose coupling of multi-physics modules into a united multi-physics system by the medium of routines listed in

Table 2. Loose coupling interface functions.

	Name	Description
A	alloc	allocation of variables x and u
I	init	initial steady-state calculation
N	next	moving forward on a time step h
S	save	saving of a state vector x
L	load	loading of the saved state vector x
U	update	update of the feedback vector u
D	free	deallocation of the variables

. For the interaction with the multi-physics system each module maintains allocation A and deallocation D routines to create and destroy a module’s internal memory space; save S and load L routines to memorize and restore a module’s state; initial steady-state I and time step marching N routines; update routine U which is used for the feedback vector update.

The feedback vectors $u_i$ are located in the external memory space shared by all the multi-physics modules, while the module’s state vectors $x_i$ are situated in the module’s internal memory spaces. This approach is employed as a multi-physics module, which usually presents a complex engineering code where the state depends on the physical variables, flags, counters, and logical descriptors needed for algorithm control. In order to avoid dealing with the code’s internal variables and mitigate the profound change of the engineering code we allow the state vectors $x_i$ to be unavailable for explicit use. The step-doubling and fixed-point algorithms can be implemented in the subspace of the feedback parameters u without an explicit operation of the corresponding state vector $x_i$.

We explain the coupling algorithm in terms of the interface functions presented in Table 2 and using the following notation: “$ABC$” denotes the successive call of the routines “A”, “B” and “C” from left to right, brackets “$\left(AB\right)$” denote a cyclic run of the routines A and B such that the sequence “$ABABAB\dots$” repeats until the convergence, the multiple dot “$A_1\dots A_n$” denotes the successive run of routines $A_1,A_2,\dots ,A_n$.

There are many methods for the solution of implicit equations such as Picard iterations, Anderson’s acceleration and Jacobi-Free Newton Krylov methods [25]. We implemented the damped Picard iteration method improved by Gauss–Seidel acceleration. The basic Picard iteration algorithm is given by

$$P_X^{\prime }=(L_1\dots L_n{X}_1\dots X_n{U}_1\dots U_n),$$

(5)

where X is the routine of initial steady-state I or time step marching N. On every iteration the multiphysics modules return to the original state by calling the subroutines $L_i$. Then routines $X_i$ update the multiphysics modules states with respect to the feedback vectors $u_i$ updated on the previous iteration by routines $U_i$. The cycle exits when the convergence for all feedback parameters is achieved.

The convergence rate of the Picard iterations can be improved by the swapping of the routines $X_i$ and $U_i$ in Equation (5) that results in a Picard iteration with the Gauss–Seidel acceleration, that is

$$P_X=(L_1\dots L_n{X}_1{U}_1\dots X_n{U}_n).$$

(6)

For the time step integration of the coupled dynamical system MPCORE uses the algorithm

$$Q_{\mathit{MPCORE}}=A_1\dots A_n{U}_1\dots U_n{S}_1\dots S_n{P}_I(S_1\dots S_n{K}_1)D_1\dots D_n,$$

(7)

where $K_1=P_N$ in case of the fixed time step method. The algorithm successively allocates the local memory space of the external modules, calls routines $U_i$ to fetch the guess values from every module and send them to the external memory space, saves the guess state of the modules using routines $S_i$, runs the Picard iteration algorithm $P_I$ to get the initial steady-state of the dynamical system, cycles time steps using $K_1$ and finally deallocates the external modules.

A time step iteration of the step-doubling method consists of two stages: the first stage performs a time step h to get less accurate solution, the second stage performs two successive time steps of size $h/2$ to get a more accurate solution. Note that $H_2$ is the routine decreasing step h into two times and $H_1$ is the routine increasing step h into two times. The one iteration of the step-doubling method is given by

$$K_2=P_N{H}_2{P}_N{S}_1\dots S_n{P}_N{H}_1$$

(8)

Time step rejection is an important feature making the multiphysics algorithm more stable to a model fail or the local truncation error (LTE) outbreak due to a large step size. In the case of such a rejection the internal memory spaces of the external models are returned to the previous state and the time step iteration repeats with a smaller step size. The drawback of the algorithms in Equations (7) and (8) is that the save routine $S_i$ is called two times in a row for each module. This results in the allocation of additional memory space for the sake of rejection capability. In order to reduce the memory requirements we use the swapped algorithm given by

$$K_3=P_N{H}_2(L_1{N}_1{N}_1{U}_1\dots L_n{N}_n{N}_n{U}_n)H_1.$$

(9)

Finally, MPCORE employs the time step integration algorithm given by (7) with $K_1$ replaced by $K_3$. Therefore, the local memory space of the external modules is saved only once per time step iteration and the rejection of the step can be achieved by the loading of external modules from the saved memory space.

2.4. Step Processing

MPCORE processes a time step according to algorithm (7) and the scheme shown in Figure 2. For a given reactor core state vector x and time step size h at the local time t, the multiphysics code advances the reactor core state from x to $x^{\prime }$ using the step-doubling and Picard iteration methods. The feedback parameters $u_k^{\prime }$ from Table 1 are evaluated for the updated solution $x^{\prime }$. The value of LTE e, for every time step is evaluated as follows

$$e=\underset{k\in \{1,\dots ,m\}}{max}{∥\frac{{\mathsf{\Delta }}_k}{A_{tol}+R_{tol}{u}_k}∥}_2,$$

(10)

where m is the number of feedback parameters, ${\mathsf{\Delta }}_k$ is the difference between two solutions obtained with the time steps h and $h/2$; $A_{tol}$ and $R_{tol}$ are the absolute and relative tolerances. For a given LTE update $e^{\prime }$, the new time step is calculated according to the step-doubling theory [13]. Note that the same norm with less relative tolerance $R_{tol}$ is used for the tracking of fixed-point iteration convergence.

The comparison of the local error, $e^{\prime }$ to 1 results in two logical branches: regular and rejected time steps. For the regular branch, $e^{\prime }\le 1$, the vector state, local time and step size are updated in order to start the next time step iteration. For the rejection branch, $e^{\prime }>1$, the solution $x^{\prime }$ is discarded and the iteration repeats with the smaller step size $h^{\prime }$ computed for the local error $e^{\prime }$. The rejection branch is used for exceptions thrown by the multi-physics code during the Picard iteration. The exception might be thrown by the Picard algorithm itself via the occurrence of either exponential divergence, periodical cycling or slow convergence rate. Moreover, the exception can be thrown by a substituent module that experience the internal divergence or failure event. For the exception of any reason the rejected iteration is repeated using the time size reduced in two times. The rejection may be repeated several times until Picard iteration converges and a small enough local error is achieved. The time step integration algorithm stops with an error if the number of rejecting cycles exceeds the limit or the time step $h^{\prime }$ becomes smaller than the minimum value.

3. Results and Analysis

MPCORE has been applied for a full-core Pressurized Water Reactor (PWR) Benchmark for Evaluation And Validation of Reactor Simulations (BEAVRS) developed by MIT Computational Reactor Physics Group [26]. The benchmark provides a detailed information about fuel assemblies, burnable absorbers, core loading patterns, in-vessel components and it enables the development and testing of detailed reactor core models for multi-physics calculation. The asymmetric control Rod Ejection Accident (REA) from Hot Zero Power (HZP) for the BEAVRS reactor core is subjected to MPCORE calculation in order to evaluate the impact of dynamic pellet-to-cladding heat transfer on design and safety of LWR.

The purpose of the coupled neutronics, thermal-hydraulics and fuel performance simulation is in more rigorous calculation of the pellet-to-cladding heat transfer to improve fuel rod temperature accuracy. Therefore, in order to evaluate the multi-physics coupling effect, the transient accident is simulated two times starting from the same initial HZP state with a constant and dynamic HTC calculated by Equation (2).

The following quantities of interest are selected for comparison: peak total core power W in GW, maximum fuel temperature $T_{fc}$ in °C, maximum cladding temperature $T_{ci}$ in °C, maximum fuel enthalpy $H_f$ in cal/g, total energy released to the coolant $Q_c$ in GJ and the number of fuel rods with the departure from nucleate boiling (DNB) event $N_{\mathrm{DNB}}$ given in% of the total number of fuel rods.

3.1. BEAVRS Core Model

The RAST-K model of the BEAVRS reactor core uses four radial nodes and 48 axial nodes per fuel assembly for the nodal diffusion calculation. The two-group cross-sections were calculated as the function of fuel temperature, coolant temperature and density, nuclear fuel burnup and boron concentration cross-section code STREAM [27] for each fuel assembly type with and without control rods.

The one-dimensional thermal-hydraulics code uses 20 axial nodes per fuel rod channel, the fuel performance code uses 10 axial nodes, 5 radial nodes in pellet and 3 radial nodes in cladding. The cladding mechanics is calculated by the one-dimensional model FRACAS-I [28], the fission gas release is calculated by FRAPFGR recommended for transient simulation, the fuel grain size is set to 10 $\mathsf{\mu }$m. The filling as-fabricated gas is helium, the cladding type is Zircaloy-4, the geometric sizes, initial enrichment and as-fabricated density of the dioxide uranium fuel depend on the location of a fuel pin with respect to the specification of BEAVRS benchmark.

3.2. HZP and HFP States

In order to evaluate the accuracy of MPCORE, the beginning of cycle (BOC) hot zero power (HZP) and hot full power (HFP) states have been verified to the benchmark data and to the high-fidelity Monte-Carlo code MCS [29]. The critical boron concentration at the BOC HZP is 937 ppm that is close to the measured value 975 ppm. The HZP state was compared to the measurement signals of radial detectors located in the central guide tubes of the fuel assemblies and showing the integral current along core elevation. The numerical signals were approximated by the arithmetic mean of the eight fuel pins neighboring a guide tube. The relative error of the numerical signals in the root-mean-square (RMS) norm accounts for about 7% while the maximum error reaches 13%. The HFP assembly-wise radial power distribution compared to MCS demonstrates the RMS relative error about 1% and 3% in the maximum. The distribution of errors in the HZP and HFP states are presented in Figure 3. Overall, we conclude the code produces an accurate distribution of power averaged over assembly; however the pin-by-pin distribution error is rather large and further improvement of the RAST-K’s reconstruction module is to be undertaken.

3.3. HZP REA Analysis

MPCORE calculates a 10-second transient at the BOC to evaluate the behavior of the BEAVRS core in the event of a prompt-critical REA. The simulation starts from the HZP core state with the total power ${10}^{-6}$ of the nominal power 3411 MWth. The control bank D at the location M-4 in benchmark coordinates has the largest worth 1377 pcm (per cent mille) or 2.0 ${\beta }_{eff}$ (effective delayed neutron fraction, 1 ${\beta }_{eff}$ = 705 pcm) as can be seen in

Table 3. Worth of control banks.

Location	Bank Type	$\mathit{\rho }$, pcm
H-8	D	49
K-8	A	175
K-6	C	319
H-2	C	360
K-2	B	589
M-4	D	1377

. This control bank is fully ejected during 100 ms at a constant rate starting from 10 ms. The transients were calculated for the relative tolerance 0.01%.

The temporal and radial distributions of the safety parameters during the transient are presented in Figure 4 and Figure 5. The reactor core is critical and all the safety parameters are constant before the control bank ejection started at 10 ms. The removal of the neutron absorber from the core produced a reactivity rise up to 2.0 ${\beta }_{eff}$. When the reactivity reached 1 ${\beta }_{eff}$ at 50 ms the reactor core turned prompt-critical and burst with a sharp pulse of neutrons. The total power exponential growth during 104 ms hit a peak of 43 nominal powers. After the peak the fuel and cladding temperatures responded immediately blowing from initial 293 °C up to 1300 °C and 520 °C respectively. The steep increase in the fuel temperature brought down the reactivity from the plateau to about 0.3 ${\beta }_{eff}$ through the Doppler effect. Therefore, the reactor core switched to delay-critical state and the total power dropped to about 1 GWth according to the amount of delay neutrons in reactor core after the power peak. The maximum fuel temperature grew gradually up to about 1550 °C during 3 s after the power peak, the pellet average enthalpy steadily decreased by means of the heat transfer from fuel to the coolant. The excessive heat leaked to the coolant and a high cladding temperature led to a boiling crisis, which was observed in the interval from 0.8 s to 2 s. The boiling crisis finished owing to the coolant circulation and the decline of the heating of the fuel.

The average gap HTC at the beginning of transient was equal to 5640 W/(${\mathrm{m}}^2$K) for both the dynamic and constant gap HTC cases. The fuel rod heating caused the thermal expansion of both a pellet and cladding materials. However, due to the thermal expansion coefficient difference, the pellet-cladding gap width diminished causing the gas HTC increment in two times with accordance to Equation (4). At the same time the radial HTC increased as well owing to a pellet temperature escalation. Overall, the change of the total gap HTC resulted in the flattening of the cladding and pellet temperatures, as shown in Figure 6. Subsequently, the outer surface cladding temperature for the dynamic gap HTC was larger compared to the constant gap HTC case increasing the likelihood of the DNB event and the number of channels experienced the coolant boiling. The fuel enthalpy, cladding temperature and fuel temperature calculated for both constant and dynamic gap HTC cases are presented in

Table 4. Peak safety parameters for constant and dynamic pellet-cladding gap HTC, $R_{tol}=0.01\%$.

Parameter	Units	Dyn	Const	$\mathsf{\Delta }$	Limits
$Wth$	$GW$	147	148	0	-
$H_f$	cal/g	90	91	1	200
$Q_c$	$MJ$	564	658	95	-
$T_{fc}$	°C	1543	1588	45	2670
$T_{ci}$	°C	531	466	−64	1649
$N_{DNB}$	%	2.9	2.5	−0.4	5

. Although the values of these parameters are below critical limits, the number of fuel rods with DNB events can violate the acceptance criteria for reactivity accident analysis if dynamic gap HTC model is in use. Note, the similar effect of the dynamic gap HTC has been reported recently for a PWR REA with a 1.3 ${\beta }_{eff}$ reactivity insertion simulated by the CTF and BISON codes [30].

The results of MPCORE simulation have been verified with the standalone RAST-K code that utilizes a nodal diffusion solver and built-in thermal-hydraulics module to produce the assembly-wise power and temperature distributions. The reactivity, total power, maximum fuel and cladding temperature functions presented in (9) demonstrate the good agreement between two codes with a maximum error about $2\%$ observed at the power and temperature peaks. The error arises primarily due to the difference in FRAPTRAN’s and RAST-K’s built-in thermal-hydraulics solvers. Note, the cross-verification of RAST-K and PARCS [31,32] codes implemented for the transient benchmarks NEACRP-PWR [33] and VVER-1000 [34] shows the difference between two codes is less 0.1%. In addition, RAST-K has been used recently for the analysis of rod ejection accident transient in a PWR core with fixed neutron source [35].

Note, although the assembly-wise distributions for both RAST-K and MPCORE codes are coincide, MPCORE produces more detailed and accurate information owing to the pin-by-pin simulation with the dynamic pellet-to-cladding heat transfer modeling.

3.4. MPCORE Performance

The convergence of the adaptive time step algorithm was analyzed using one of the BEAVRS’ fuel assemblies. The fuel assembly under consideration acquires enrichment 1.6%, coolant flow rate $G=84.1$ kg/s and initial total power 17.7 Wth; the assembly is modeled with the reflective boundary conditions in the radial directions and zero-current boundary conditions at the assembly top and bottom where reflector is placed. The power pulse was initiated by the reduction of boron concentration from the critical 905 ppm to sub-critical 825 ppm during 10 ms, which is equivalent to the insertion of reactivity 1349 pcm. According to RAST-K calculation, conducted for the constant gap HTC 5640 W/${\mathrm{m}}^2$K and time step $h=0.1$ ms, the power peak of the transient occurs at the time 191.9 ms and accounts for 2.220 GWth; this value was set as a reference for further convergence study.

The fuel assembly was simulated by MPCORE with the constant and dynamic gap HTC models using various relative tolerances $R_{tol}$ taken between ${10}^{-3}$ and 1%, while the relative tolerance used for Picard iteration is 100 times less. The graphics of the MPCORE’s execution times and power peak errors are presented in Figure 7. The constant and dynamic models show exponential behavior on the relative tolerance; the dynamic case consumes more computation time owing to slower convergence of Picard iterations. The relative number of rejected time steps varies between 3% and 9%. The power peak value converges to approximately the same number for both dynamic and constant gap HTC indicating the initial stage of a REA transient is almost insensitive to the pellet-to-cladding heat transfer simulation.

The analysis of the multi-physics code performance for the full-core case has been done for $R_{tol}=0.01\%$, the results are presented in Figure 8. Overall, the simulation completed 363 regular and 55 rejected time steps, 3340 internal iterations and consumed about 42 h of execution time using one computing processor. The graphic of time step size history shows the minimum step size 1 ms was at the power peak region; afterwards the step size gradually increased accounting for about 0.5 s at the end of the simulation interval. The number of function calls per time step iteration was equal to 6 during power pulse; then it increased up to 10 due in response to a thermal-hydraulics model. Some of time steps required over 30 iterations to be converged, which indicates the slow convergence of the Picard iteration applied for the chosen step size. Note, that each iteration includes one large and two small time steps performed according to the step-doubling approach. The number of function calls remains the same for all substituent modules.

The most consuming part of the multi-physics simulation is fuel performance that accounts over 50% of the total execution time. The thermal-hydraulics module requires about 15% during the calculation of power pulse and 30% at the end of the simulation interval. The neutronics calculations consume less than 15% over all time. Note, that MPCORE implements integration of neutron kinetics, thermal-hydraulics and fuel performance modules according to algorithm Equation (7) with the data exchange through a shared memory space. Therefore, the memory-handling routines including $S_i,L_i$, $U_i$ and damping data routine can be slow as the rank of feedback vectors is large; in this study the overall consumption of the memory-related routine is about $8\%$.

4. Conclusions

The nodal diffusion code RAST-K, fuel performance analysis codes FRAPTRAN, and one-dimensional coolant thermal-hydraulics code CTH1D have been coupled together by means of code MPCORE for the high-resolution multi-physics analysis of transient PWR cores. The new multi-physics system implements the loose coupling of constituent modules enhanced by an adaptive time step algorithm for solution accuracy control and by the restart capability for robustness improvement. Note, each module can be replaced with little efforts on a more suitable or accurate one owing to the modular structure and transparent coupling interface supported by the new coupling system. The present configuration of MPCORE comprised of the two-group nodal diffusion, one-dimensional thermal-hydraulics and thermal-mechanics models is dedicated to the core design optimization and uncertainty quantification where the computational performance is highly demanded for reactor core multiple runs.

The coupled system of multi-physics equations is integrated in time using the step-doubling method to achieve automatic step size and the damped Picard iterations to resolve the non-linear system of equations associated with a given time step. The convergence of the Picard method has been improved by the Gauss–Seidel acceleration algorithm in which a multi-physics module is fulfilled with respect to the feedback information newly updated by a preceding module in the execution chain. Furthermore, the random access memory requirement has been reduced by swapping of the time step and Picard iteration loops in the total coupling algorithm to avoid the allocation of an extra backup memory space needed for the time step restart option.

The new multi-physics system has been applied for the simulation of a hypothetical severe transient accident initiated by ejection of control rods in a PWR core. The BEAVRS benchmark core has been chosen for the test as it provides comprehensive information about fuel pins and assemblies, and enables to develop a detailed multi-physics reactor core model. The simulation was performed in order to assess the accuracy, convergence and performance of MPCORE, as well as in the interest of the effect of dynamic pellet-to-cladding heat transfer modeling on reactor core safety parameters.

Overall, it was found that, while the peaking values of the total core power and the fuel enthalpy are practically indifferent to the choice of constant or dynamic heat transfer models, the impact on the fuel and cladding temperatures is much more considerable. In detail, the increase in the heat transfer coefficient, observed in the dynamic model due to the pellet-cladding gap closing, facilitates the transfer of heat from a fuel pellet to cladding. The rise in the cladding outer surface temperature makes DNB events happen more likely for the larger number of fuel pins. Therefore, the simulation performed using the dynamic gap HTC model can potentially breach the acceptance criteria for reactivity accident analysis, although the fixed gap HTC model produces satisfied results. Note, that the fixed model can be made more conservative using a larger value of HTC than the value fixed at the beginning of transient. However, the conservative HTC value is still unknown and requires either the dynamic heat transfer modeling or an expert’s judgment.

The future work will be focused on the improvement of the inner iterations convergence, and on the development of the cladding-to-fluid heat transfer performance. The full-core simulation reveals the slow convergence of Picard iterations with a fixed damping factor resulting in increase of the rejected time steps number and decrease in a time step size. On the other hand, as was reported in [25], the considerable improvement in convergence rate can be achieved using Anderson’s acceleration method that requires a minor modification of the existing iteration scheme. Notwithstanding the relatively slow convergence, the accomplished full-core multi-physics analysis has demonstrated significance of the pellet-to-cladding heat transfer modeling for a severe nuclear accident in a PWR core. However, further investigation of the influence of the model parameters and resolutions on the target reactor core safety characteristics is of vital importance.