BEDOK: An In-House Numerical Reactor Simulator Effort in Singapore

1  Introduction

Nuclear power produces a large amount of electricity with a relatively small environmental footprint, providing a stable and reliable power supply while mitigating climate change. Thus it is seen as essential for meeting the growing energy demands of modern societies. Modern nuclear reactors are designed with multiple layers of redundant safety systems to prevent accidents and manage potential risks. However, nuclear reactor experiments are complex, long-term, costly and highly sensitive; it is impractical to verify the reliability and safety of nuclear reactor design solely through experiments. High-performance computing (HPC) technology that supports high-fidelity numerical simulation models has become an indispensable way to study and develop nuclear reactors. At present, developed countries around the world, such as the United States, China and the EU as a whole have carried out extensive research in the field of numerical reactors. For example, the United States has launched three large-scale numerical reactor research projects: the Consortium for Advanced Simulation of Light Water Reactors (CASL) [1], Nuclear Energy Advanced Modeling and Simulation (NEAMS) [2], and Center for Exascale Simulation of Advanced Reactors (CESAR). These projects are aimed at the safety analysis of existing light water reactors (LWR) and next-generation reactors as well as the development of E-class supercomputers to carry out advanced numerical reactor research. Similarly, China has also developed an exascale simulator, China Virtual Reactor (CVR) [3,4], that aims to achieve full-core high-fidelity simulation. In the EU, the CORTEX (CORe monitoring Techniques and EXperimental validation and demonstration) and the McSAFER (High-Performance Advanced Methods and Experimental Investigations for the Safety Evaluation of Generic Small Modular Reactors) projects develop neutronic, thermal-hydraulic, and thermo-mechanic simulators [5]. On the other hand, prospective adopters of nuclear energy, such as the UAE, have also started development of their own numerical simulator projects [6]. In addition to national initiatives, private operators of nuclear power plants have also been making significant investments in their own numerical simulator platforms [7].

Numerical simulation research is therefore crucial for countries that are considering the adaptation of nuclear power. It is in this context that this work aims to develop the framework for a long-term numerical simulator development project under the title Broad-scope Environment for Dual-phase advanced reactor Operation simulation Kit (BEDOK). A compilation of coupled numerical codes is to be developed in-house, encompassing each of the major topics in numerical reactor simulation, including cross-section generation, core neutronics, thermal hydraulics, fuel performance calculation and multiphysics coupling. The development of in-house simulation capabilities provides Singapore with independent methods of reactor safety assessment and analysis. BEDOK follows previous in-house efforts in monte-carlo neutronics [8,9], as well as in a molten-salt GUI simulator [10].

The current work presented here represents the first step in BEDOK, in which we start from lower fidelity methods initially, which encompasses a 3D steady state neutron flux solver based on the semi-analytical nodal diffusion method. Higher fidelity methods will be developed in the future as replacements. The neutronics solver is coupled with a 1D thermal hydraulic (T-H) solver with cross-section feedback. The development of BEDOK will prioritise modularity as well as the framework to readily allow for ongoing incorporation of additional capabilities as well as switching to higher fidelity models. This is essential due to rapid developments in the field of numerical reactor simulation. Each current capability in BEDOK has also been verified against benchmark problems. Neutron flux output in this work has been benchmarked against the standard IAEA 2D and 3D PWR benchmark [11], while neutronics coupling with T-H parameters in this work is verified via the NEACRP 3D LWR Core Transient Benchmark [12,13]. The 1D thermal hydraulic drift-flux model is verified against the OECD/NRC benchmark based on NUPEC PWR Sub-channel and Bundle Tests (PSBT) [14,15] Test Series 1.

The rest of the paper is organised as follows. In Section 2, the theory behind each solver is introduced, as well as the benchmark problems used in the verification and validation of the constituent models in BEDOK. In Section 3, results of each benchmark case are reported, with comparison to results obtained by other benchmark participants as well as other simulators that also have published results for the aforementioned benchmark cases. Section 4 summarizes the current progress of BEDOK as well as outlining future improvements.

2  Models and Framework

2.1 Neutronics Model

The neutron transport equation can be integrated over Ω and energy groups g to obtain the diffusion approximation as follows

1vg∂∂tϕg(x,t)=−∇⋅Jg(x,t)−σrg(x,t)ϕg(x,t)+Qg(x,t)(1)

where the energy and angle integrated source term is

Qg(x,t)=1λ∑g′χgνσfg′(x,t)ϕg′(x,t)+∑g′≠gσsgg′(x,t)ϕg′(x,t)(2)

where λ is an eigenvalue, χ is the neutron generation fraction, ν is the neutron multiplicity. σr, σf and σs refer to the transport, fission and scattering cross sections, respectively.

J(x,t)=∫Ωψ(x,Ω,t) dΩ(3)

is the neutron current. The energy group index g is omitted for convenience. Then the approximation

J(x,t)≈−D(x)∇ϕ(x,t)(4)

for diffusion coefficient D comes from Fick’s law, with the scalar flux ϕ defined as

ϕ(x,t)=∫ψ(x,Ω,t) dΩ(5)

The neutron diffusion can discretised over each equally spaced node in the problem geometry to obtain then be integrated over two of the three dimensions to obtain 1D diffusion equations (shown here for x-direction) coupled via additional transverse leakage source terms:

Lx(x,t)=1ΔyΔz∫∫Jy(x,t)+Jz(x,t) dydz(6)

where Δ refer to the discretised cell lengths and Jy and Jz refer to the respective directional currents. The basis of the nodal diffusion method then comes when solving for ϕ at each node. For an in depth classical review of nodal methods in neutron transport, the reader is referred to work by Lawrence [16]. In nodal expansion methods, 1D fluxes are further expanded separately for each node. In the case of a fourth-order semi-analytical nodal method (SANM) [17], the expansion functions are:

f0(x)=1f1(x)=xf2(x)=12(3x2+1)f3(x)=sinh⁡(βgx)−ℳg1(sinh)f1(x)sinh⁡(βg)−ℳg1(sinh)f4(x)=cosh⁡(βgx)−ℳg0(cosh)f0(x)−ℳg2(cosh)f2(x)cosh⁡(βg)−ℳg0(cosh)−2ℳg2(cosh)(7)

with

βg=ΣtgDgΔx2(8)

and

ℳgk(F)=22k+1∫−11F(βgx)fk(x) dx(9)

The resultant group discretised scalar flux is given by

ϕg(x)=ϕ¯g+∑n=14agnfn(2xΔx)(10)

where ϕ¯g is the average scalar flux over the node and the an are derived separately for each node [18,19]. The relations for an derived in BEDOK follows work by Tang [19], which is related to the boundary flux currents, leakage moments and buckling terms as

ag1=ag11(Jg++Jg−)+ag12(∑g≠g′Bg′gag′1)−Lg1ag2=ag21(Jg++Jg−)+ag22(∑g≠g′Bg′gag′2)−Lg2ag3=ag31(Jg++Jg−)+ag32(∑g≠g′Bg′gag′1)−Lg1ag4=ag41(Jg++Jg−)+ag42(∑g≠g′Bg′gag′2)−Lg2(11)

The buckling terms are given by

Bgg=ΣtgandBgg′=−Σgg′−χgkeffνΣfg(12)

where χ is the neutron generation spectrum and ν is the neutrons generated per fission. The boundary currents are

Jg+=−2DgΔx(ag1+3ag2+ℋgag3+𝒢gag4)(13)

and

Jg−=−2DgΔx(ag1−3ag2+ℋgag3−𝒢gag4)(14)

where

ℋg=βgcosh⁡(βg)−ℳg1(cosh)sinh⁡(βg)−ℳg1(sinh)(15)

and

𝒢g=βgsinh⁡(βg)−3ℳg2(cosh)cosh⁡(βg)−ℳg0(cosh)−ℳg2(cosh)(16)

Lgn are the leakage moments obtained from a second order expansion of the transverse leakage source term with

Lg(x)=Lg0+∑n=12Lgnfn(2xΔx)(17)

The agnn are components of a ann vector over the group indices, derived by

a11=−Δx4D(1+ℋ𝒜B)a12=ℋ4D(1+ℋ𝒜B)⋅𝒜a21=−Δx4D(3+𝒢𝒞B)a22=𝒢D(3+𝒢𝒞B)⋅𝒞a31=−Δx4D(1+ℋ𝒜B)⋅𝒜Ba32=−14D(1+ℋ𝒜B)⋅𝒜a41=−Δx4D(3+𝒢𝒞B)⋅𝒞Ba42=−3D(3+𝒢𝒞B)⋅𝒞(18)

The remaining functions are defined as

𝒜g=sinh⁡(βg)−ℳg1(cosh)Σtgℳg1(sinh)(19)

and

𝒞g=cosh⁡(βg)−ℳg0(cosh)−ℳg2(cosh)Σtgℳg2(sinh)(20)

With the functions determined, Eqs. (13) and (14) can be substituted back into Eq. (11) to solve for the coefficients an via matrix inversion.

2.2 Thermal Hydraulic Model

Thermal hydraulic feedback to the neutronic system is achieved via cross section variations with respect to boron density, moderator temperature, moderator density, and fuel temperature. For property λ the cross sections vary as

σ=σo+∂σ∂λ|o(λ−λo)(21)

about reference point “o”, with the exception of fuel doppler temperature which varies the cross sections with

σ=σo+∂σ∂λ|o(λ−λo)(22)

Coolant temperatures and densities are obtained from a simple constant pressure 1D model where coolant enthalpy is integrated from the inlet as

H(z)=∫q′′(z)GmAf dz(23)

where q′′ is the heat flux per unit area, Gm is the coolant mass flow rate and Af is the flow area. Other thermodynamic properties such as temperature and density are subsequently recovered via look-up table generated from IAPWS-IF97 data [20]. The grid size of the lookup table required for accurate results has been discussed by Zou et al. [21] and is adapted in this work. Temperatures in the fuel rod are obtained via 1D heat equation assuming even power generation within the fueled rod cross section.

In addition to the single phase coolant model, and alternative two-phase model based on the four-equation drift flux model is available based on work by Zou et al. [22]. The constituent equations of the drift flux model are firstly the mixture (subscript m) continuity equation

∂ρm∂t+∂(ρmvm)∂z=0(24)

The second is the continuity equation for the dispersed phase (i.e., the gas phase g),

∂(αgρg)∂t+∂(αgρgvm)∂z+∂∂z(αgρgρlρmV¯g)=Γg(25)

where subscript l refers to the liquid phase, Γg and V¯g is the generation rate and mean drift velocity of phase g. Next is the mixture momentum equation

∂vm∂t+vm∂vm∂z=−1ρm∂p∂z−g→⋅v→m|v→m|−fm2DHvm|vm|−1ρm∂∂z(αgρgρlαlρmV¯g2)−1ρm∂∂z(τzz+τ~zz+∑k⟨αkρkvk(vk−⟨αkvk⟩⟨αk⟩)⟩)(26)

where f is the friction factor, p is the pressure and DH is the hydraulic diameter, τ is the normal viscous stress tensor while τ~ is the normal turbulent stress tensor and ⟨ ⟩ refers to a flow area averaged value. Finally, there is the mixture enthalpy-energy equation

∂(ρmEm)∂t+∂(ρmhmvm)∂z+∂∂z(αgρgρlρmV¯g(⟨αghg⟩⟨αg⟩−⟨αlhl⟩⟨αl⟩))=qw′′′aw+vm∂p∂z+αg(ρl−ρg)ρmV¯g∂p∂z(27)

where E is the internal energy, h is the enthalpy, qw′′′ is the wall heat flux per unit volume and aw is the heated surface area. The variables to be solved for are p, αg, vm and T. Terms such as V¯g and Γg are implicitly defined from the input variables with a functional form that can depend on the type of flow regime, and therefore would constitute separate piecewise models in such cases. The flow regimes considered in this work are bubbly flow, slug flow, annular mist flow and a transition region between the slug and annular flow regime. The flow regime regions are separated by cutoff αg values as shown in Fig. 1, which are implemented with the same criteria as in RELAP5/MOD3 [23]. The functional forms of V¯g and Γg in each regime are referenced from TRACE [24] and RELAP5/MOD3 manuals [23], as well as work by Zou et al. [22].

Figure 1: Flow regime map used in this work, which is dependent on gas fraction αg.

In this work the void generation rate Γg models are referenced from the RELAP5/MOD3 manuals [23]. The models are indeed piecewise according to Fig. 1 and thus not shown here due to the complexity.

While piecewise models for V¯g are available [25–29], the ERPI drift-flux closure correlations [30], which consist of only one functional model, were validated against the full relevant range of flow regimes and also subsequently modified and adapted in RELAP5 [23]. The EPRI model for V¯g is given by

V¯g=ρmVj+ρm(C0−1)vmρm−(C0−1)αg(ρl−ρg)(28)

with

C0=L0K0+(1+K0)αgr0(29)

being the distribution parameter and

Vj=1.41((ρl−ρg)gσρl2)0.25C2C3C4C5(30)

σ is the surface tension for which this work uses the model from Vargaftik et al. [31] with

σ=235.8(1−0.625(647.15−Tl)647.15)(647.15−Tl647.15)1.256(31)

L0, K0, r0 as well as C2 to C5 are empirically fitted functions of the primary drift flux variables p, αg, vm and T which follow the description in the RELAP5 manuals [23].

The friction factor f in Eq. (26) is determined from wall friction pressure drop with

(∂p∂z)m=ϕg(∂p∂z)g=ϕl(∂p∂z)l(32)

where ϕ refers to the Darcy-Weisbach friction multipliers. This then relates to the friction factor through

fmvm|vm|=1ρm(λlρlαl2vl2+Cfλlρlαl2vl2λgρgαg2vg2+λgρgαg2vg2)(33)

where

Cf=−2+(28−0.3Gm)exp⁡((log10⁡Λ+2.5)22.4−0.0001Gm)(34)

and

Λ=ρgρl(μlμg)0.25(35)

μ is the dynamic viscosity, Gm is the mass flow rate and λ is the Darcy friction factor, for which this work employs the model developed by Cheng [32].

The spatial discretisation for the drift flux problem follows a first order donor-cell scheme. Within this scheme, derivatives are taken via finite differencing of ‘upstream’ nodes. For instance considering the derivative for an arbitrary collection of terms F

∂Fi∂z=1Δz(Fi+1/2∗−Fi−1/2∗)(36)

then if vm>0, velocity terms will be evaluated as per normal while other terms such as α and ρ in F∗ are evaluated upstream at i instead of i+1/2 and F1/2 will be the inlet values. For example, in the mixture continuity equation

∂(ρmvm)∂z|i=1Δz,i(ρm,ivm,i+1/2−ρm,i−1vm,i−1/2)(37)

Variables evaluated at i+1/2 are taken as the average value of the neighbouring nodes, for example

ρm,i+1/2=12(ρm,i+ρm,i+1)(38)

Boiling Crisis

In LWRs, the critical heat flux (CHF) is a important safety topic concerning the heat flux threshold at which surface boiling ceases to be an effective form of transferring heat into the coolant. In this work, we analyse the CHF of PWR systems by employing the Westinghouse-3 (W-3) correlations [33], given by

qCHF′′=K1K2K3K4(39)

where

K1=(2.022−0.06238p)+(0.1722−0.01427p)exp⁡((18.177−0.5987p)χe)K2=(0.1484−1.596χe+0.1729χe|χe|)2.326Gm+3271K3=(1.157−0.869χe)(0.2664+0.8357exp⁡(−124.1DH/100))K4=(0.8258+0.0003413(hl,sat−hin))(40)

where hl,sat is the saturated liquid enthalpy and hin is the inlet enthalpy for the node segment. χe is the equilibrium thermodynamic quality given by

χe=h−hlhlg(41)

where hlg is the latent heat of vaporisation. Following which, the relevant figure is the departure from nucleate boiling ratio (DNBR), which is the ratio between predicted CHF value and the actual heat flux from the fuel rod.

In BWRs, the relevant parameter used to analyse the safety margins with respect to boiling crisis is the critical power ratio (CPR). This is typically calculated from the critical quality χcrit which is where dry-out is predicted to occur. In simulators the thermal power can then be increased till χcrit is reached in the core. This increased thermal power ratio is the CPR. This work employs the M3-CISE4 correlation [34] for χcrit which for a given channel is

χcrit=(DHDq)23K5LqK6+Lq(χ¯χmax)0.3(42)

where

K5=[1+0.0001481(1−PPcrit)−3Gm]−1forGm⩽3375(1−PPcrit)3kg/s/m2(43)

K5=(1−PPcrit)(1000Gm)13forGm>3375(1−PPcrit)3kg/s/m2(44)

and

K6=0.279(PcritP−1)0.4GmDH1.4(45)

Pcrit is the critical pressure, Dq is the heated parameter and Lq is the heated length.

2.3 Code Structure and Benchmarking

The overall structure and implementation of BEDOK is shown in Fig. 2, which illustrates the solver modules and data structures passed between them. Input data is passed to a neutronics solver of choice, which outputs the power profile to the T-H solver. The result is then used to update the cross sections, which then feeds back to the neutronics solver. The process is repeated till convergence is reached. The modular nature of this work means that different models will be developed for each solver and should pass the same data types to the next step in the process. Within each solver type, it is also possible to use the lower fidelity but computationally faster models as acceleration for higher fidelity models.

Figure 2: Flowchart for the overall algorithm in this work, showing each separate module and data type passed between them.

Further, iterative systems are accelerated with the use of an extended Anderson scheme [35] based on the original Anderson acceleration [36]. A range of acceleration schemes for user choice is also intended to be implemented in BEDOK. The extended Anderson scheme is tailored for fixed-point iteration problems but can be easily modified to solve minimization problems as well. For solution iterate Xi at iteration i, a sequence of of coefficients ak is calculated, for a scheme using history length L such that the next iterate is given by

Xi+1=Xi+βRi−∑m=i−Li−1am(Xm+1+βRm+1−Xm−βRm)(46)

where R is the residual and β is the relaxation factor. A optimised solution update is thus obtained using information from the past L iterations. The Anderson scheme is also known to considerably reduce coupling iterations between the neutronic and T-H systems by smoothing local oscillatory behaviors [37], and is able to handle problems where ill-conditioned matrices can cause convergence issues in Newton methods [36]. This is especially relevant here as the non-linear drift flux equations are solved using the Jacobian-free Newton–Krylov (JFNK) method [38].

This work is written in MATLAB and simulations are executed in parallel using 10 cores on the Intel Xeon W-2255 CPU @ 3.70 GHz processor. Neutron flux output in this work is first benchmarked against the standard IAEA 2D and 3D PWR benchmark [11], with specifications also listed in Ref. [39]. Neutronics coupling with T-H parameters in this work is verified via the NEACRP 3-D LWR Core Transient Benchmark [12,13], in particular, case A2 of the PWR benchmarks. A quadrant geometry is used with reflective boundaries to simulate the full core. Further void fraction results of the drift-flux model are verified against the OECD/NRC benchmark based on NUPEC PWR Sub-channel and Bundle Tests (PSBT) [14,15] Test Series 1. The experiment behind the Test series 1 benchmark involves a high-pressure and high-temperature test loop simulating various sub-channel types found in a PWR assembly. The heated section of the loop also contains a void measurement section using gamma-ray transmission method.

3  Results and Analysis

3.1 PWR Benchmarking

The steady state IAEA PWR keff are benchmarked to results obtained from PARCS [40], a nodal diffusion solver developed under the US NRC. The results are presented in Table 1. keff obtained via pure diffusion solver are within 0.01% difference with PARCS while the nodal diffusion method yields a difference of 0.001% for both 2D and 3D cases, a significant improvement as expected of a higher fidelity method. The core power distribution of the IAEA 3D benchmark for each x-y cell is also presented in Fig. 3, compared against available data from the open nuclear reactor simulator KOMODO [6], which present core power densities within 1% difference. The 3D power profile of the IAEA 3D benchmark case is also generated and shown in Fig. 4, where it can be seen that due to a lack of a T-H system feedback, the power is profile is concentrated in the center of the core with regions of lower power corresponding to presence of control rods.

Figure 3: Core power distribution for IAEA 3D benchmark in comparison to KOMODO results [6]. The upper values in each cell are the KOMODO values and the lower ones are from this work.

Figure 4: 3D core power profile for the IAEA 3D benchmark [11] generated in this work.

With the accuracy of the nodal diffusion solver in BEDOK verified, the next step is to couple the neutronic system to a T-H system that is driven by 1D solvers for both fuel and coolant temperature. For the T-H coupled case, the benchmark of choice is case A2 of the NEACRP PWR benchmarks [11]. The keff results are presented in Table 2 with reference values. The critical boron concentration with no T-H feedback assumes fuel and coolant temperatures at reference values. It is clear that the addition of coupling increases the error of reference keff value to 1.7% in this work. As the T-H system employed in this work consists of lower fidelity 1D models, is expected that the error can be reduced further with more sophisticated T-H models. It is again useful to note here that the modular nature BEDOK is intended to allow efficient switching of individual models, such as is expected for the T-H model here, to be crucial for the work in progress. It can then be seen from Table 3 that slight improvement is obtained from using a two-phase model, though difference is small for as void fractions are expected to be very small in PWR systems. The 3D power profile of case A2 is also illustrated in Fig. 5, where the difference with the simpler IAEA 3D case shown in Fig. 4 can be seen with a lower axial power peak. This derives from the negative reactivity coefficient of the core driving down power as the coolant heats up. The central control rod region is also distinguishable in the 3D power profile.

Figure 5: 3D core power profile for the A2 benchmark case in the NEACRP 3-D LWR Core Transient Benchmark [12,13] generated in this work.

When the analysis is performed via the W-3 CHF correlation, the minimum acceptable DNBR is 1.3 [41]. For the PWR benchmark base A2, the predicted DNBR is plot in Fig. 6 alongside the calculated CHF, the simulation heat flux, as well as the minimum DNBR. The minimum calculated DNBR for this case is 2.4, well above the minimum allowed DNBR, as expected of a PWR in normal operation.

Figure 6: Predicted departure from nucleate boiling ratio (DNBR) for the A2 benchmark case in the NEACRP 3-D LWR Core Transient Benchmark [12,13]. The critical heat flux, calculated heat flux ad minimum acceptable DNBR of 1.3 is also displayed.

3.2 Two-Phase Flow Benchmarking

The OECD/NRC PBST benchmark [14,15] includes participants using computational fluid dynamic codes of various complexity and fidelity such as RELAP5 [23], TRACE [24] NEPTUNE [43], CATHARE-3 [44] and FLICA-4 [45]. The four equation drift-flux model developed in this work is benchmark to the Test series 1 of the OECD/NRC PBST benchmark [14], with the obtained void fraction values presented in Table 4. Figs. 7–11 are graphs adapted from the NEACRP benchmark results [13], which gives a comparison of how the αg profile generated in this work compares to results obtained by the various benchmark participants. For the most part, it can be seen that this work performs fairly well in comparison, though it is noted that Run 1.2211 that features a low experimental void fraction value, where this work predicts a void fraction of 0.2 in comparison to the experimental value of 0.038 at axial level 1.4 m. It is also noted from Fig. 7 that while may other benchmark participants also significantly over predict the void fraction value, even including commercial 3D computational fluid dynamics codes, this work still predicts the highest αg. This signifies some room for future development on the low void fraction regime. Regardless, this work performs on par with the benchmark participants for the remaining test cases, as well as obtaining the most accurate result in Run 1.2237 comparing to other work in Fig. 9. Further Run 1.2221 also features a low void fraction test case for which the current model is able to accurately predict. Runs 1.5221 and 1.6221 feature similar low void fraction cases where the current model over predicts the void fraction.

Figure 7: Axial void fraction profile for Run 1.2211 of the OECD/NRC PBST benchmark [14] in comparison with experiment and benchmark participants, graph adapted from the result publications [15].

Figure 8: Axial void fraction profile for Run 1.2223 of the OECD/NRC PBST benchmark [14] in comparison with experiment and benchmark participants, graph adapted from the result publications [15].

Figure 9: Axial void fraction profile for Run 1.2237 of the OECD/NRC PBST benchmark [14] in comparison with experiment and benchmark participants, graph adapted from the result publications [15].

Figure 10: Axial void fraction profile for Run 1.4325 of the OECD/NRC PBST benchmark [14] in comparison with experiment and benchmark participants, graph adapted from the result publications [15].

Figure 11: Axial void fraction profile for Run 1.4326 of the OECD/NRC PBST benchmark [14] in comparison with experiment and benchmark participants, graph adapted from the result publications [15].

3.3 Drift-Flux Convergence

Due to nonlinear closure models in the four equation drift flux model, difficulties were expected in the numerical convergence of the two-phase flow solutions. In this aspect, various techniques to achieve better numerical stability such as second order differencing and staggered grid discretisation has been studied [21,22]. In general, various methods are available for the acceleration and convergence improvements of the solution of non-linear equations, especially those derived from physical systems [38,46–50]. Not all methods have been applied specifically to the solution of drift-flux equations. The variety of approaches available, such as reduced order methods for example [49,50], mean that the actual effectiveness of these techniques in the context of the drift-flux equations would require extensive effort and thus be more appropriate as a future publication. A more concise effort is mode in this work where the effects of a few simpler stabilisation and acceleration techniques are studied in this work to demonstrate the possibilities available for improving the computational efficiency of the drift-flux equations. The first is varying the weights of the drift flux equations and the second is the utilisation of the extended Anderson acceleration scheme [35]. Preconditioning of the JFNK solver is not done as implementing the finite difference preconditioner in MATLAB was unsuccessful.

The varying of equation weight is implemented as follows, the four drift-flux equations in BEDOK are implemented with base units g/cm3/s for Eqs. (24) and (25), cm/s2 for Eq. (26) and g/cm/s3 for Eq. (27). Eqs. (24)–(27) are then each multiplied by a weight factor wi for i=1,2,3,4. So for instance for Eq. (24), the associated residue ℛ1 is

∂ρm∂t+∂(ρmvm)∂z=ℛ1(47)

which after including the weight factor becomes

w1[∂ρm∂t+∂(ρmvm)∂z]=ℛ1(48)

To ensure that faster convergence is not being achieved due to the weight factors relaxing the convergence criteria, the condition wi⩾1 is imposed and convergence will be judged based on the sum of the L2-norm of the equation residues. A non-exhaustive study is done on Run 1.2237 of the OECD/NRC PBST benchmark [14] with results reported in Table 5. It can be seen that the base case of w1=w2=w3=w4=1 actually does not converge, but increasing the weights appropriately does allow convergence of the equations even when the convergence criteria is strictly tighter. The case of w1=10, w2=5, w3=1, w4=100 performed the best here and is thus chosen for BEDOK, but a more exhaustive study could perhaps yield a better weight distribution, and the trends reflected in Run 1.2237 may not necessarily hold in general.

The effectiveness of the extended Anderson acceleration scheme [35] is studied in the context of a subset of cases in the test series 1 of the OECD/NRC PBST benchmark [14], with the convergence behaviours reported in Figs. 12–16. In all cases it can be seen that the extended Anderson acceleration scheme does provide faster convergence in comparison to cases where the scheme is not used. It is also noted that the accelerated cases tend to have stepwise convergence behaviour as opposed to a slower but steady convergence rate of the unaccelerated case.

Figure 12: Convergence behaviour for Run 1.2211 of the OECD/NRC PBST benchmark [14] with and without the extended Anderson scheme [35]. Anderson-n indicates the history length used.

Figure 13: Convergence behaviour for Run 1.2223 of the OECD/NRC PBST benchmark [14] with and without the extended Anderson scheme [35]. Anderson-n indicates the history length used.

Figure 14: Convergence behaviour for Run 1.2237 of the OECD/NRC PBST benchmark [14] with and without the extended Anderson scheme [35]. Anderson-n indicates the history length used.

Figure 15: Convergence behaviour for Run 1.4325 of the OECD/NRC PBST benchmark [14] with and without the extended Anderson scheme [35]. Anderson-n indicates the history length used.

Figure 16: Convergence behaviour for Run 1.4326 of the OECD/NRC PBST benchmark [14] with and without the extended Anderson scheme [35]. Anderson-n indicates the history length used.

3.4 BWR Benchmarking

With a two-phase T-H model and a neutronics solver, the next step is to verify the coupling between the two via the NEACRP BWR case D1 benchmark. The steady state keff value obtained is indicated in Table 6, where this work obtains a keff value within 0.2% of the benchmark participants, well within the standard deviation. The core power distribution for this BWR benchmark is also presented in Fig. 17. Following, the critical power ratio for this BWR benchmark case is also analysed for which in this work the critical void fraction predicted by the MS-CISE4 correlation [34] is presented in Fig. 18 in relation to the hottest channel, which shows a safety margin expected in steady state operation.

Figure 17: 3D core power profile for the D1 benchmark case in the NEACRP 3-D LWR Core Transient Benchmark [12,13] generated in this work.

Figure 18: Predicted critical void fraction, based on M3-CISE4 critical equilibrium correlation [34], in the D1 benchmark case in the NEACRP 3-D LWR Core Transient Benchmark [12,13], in comparison to the hottest channel.

4  Conclusion

The development of Singapore’s in-house nuclear reactor simulator code package, BEDOK, marks a first step in Singapore’s commitment in advancing local expertise in the development and understanding of nuclear reactor simulation codes. This simulation tool can currently be applied to analysis in steady-state LWR problems, with further development ongoing. BEDOK currently supports a nodal diffusion solver for neutronics that employs a fourth-order expansion on the flux. Whereas for the T-H module, a single-phase 1D model is included for analysis in PWRs while a 1D drift-flux model was also developed for two-phase flow analysis. Capabilities of BEDOK are demonstrated in the verification of IAEA [11] and NEACRP [12,13] benchmark cases for generic LWR reactors, as well as OECD PWR sub-channel and bundle tests for two-phase flow verification [14,15]. The benchmark verification provides confidence in the core neutronics driver, the T-H models, as well as the coupling regime for these two systems in time-independent problems. Some difficulties were encountered in the convergence of the drift flux solutions, as discussed with some novel insights. The implementation of further stabilisation techniques will be included in the ongoing development of BEDOK. Parameters critical to reactor safety margins, such as CHF and DNBR have also been calculated in the benchmark cases. Future work will include development of time-dependent capabilities, optimisation of algorithms for efficiency, and further refinement of accuracy. As part of the BEDOK project, cross-section generation capabilities are also to be developed to enable the analysis of new-generation small modular reactors, which are the reactor type of key interest in the Southeast Asian region.

Acknowledgement: The authors thank the Singapore Nuclear Research and Safety Institute, National University of Singapore for computational resources provided.

Funding Statement: This work was supported by the National Research Foundation Singapore (A-8002968-00-00).

Author Contributions: The authors confirm contribution to the paper as follows—Yan Ren Than: writing (original draft preparation, reviewing and editing), conceptualisation, code writing and validation, investigation, formal analysis. Sicong Xiao: writing (reviewing and editing), conceptualisation, supervision, formal analysis. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The data that supports the findings of this study are available from the corresponding author upon reasonable request.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Turinsky PJ, Cothe DB. Modeling and simulation challenges pursued by the consortium for advanced simulation of light water reactors (CASL). J Comput Phys. 2016;326(8):544–68. doi:10.1016/j.jcp.2016.02.043. [Google Scholar] [CrossRef]

2. Sofu T, Thomas JW. US DOE NEAMS program and SHARP multi-physics toolkit for high-fidelity SFR core design and analysis. In: Proceedings of the International Conference on Fast Reactors and Related Fuel Cycles: Next Generation Nuclear Systems for Sustainable Development (FR17); 2017 Jun 26–29; Yekaterinburg, Russia. [Google Scholar]

3. Lu X, Li Y, Chen D, Chu G, Wang A. Challenges of high-fidelity virtual reactor for exascale computing and research progress of China virtual reactor. Nucl Eng Des. 2023;413(3):112566. doi:10.1016/j.nucengdes.2023.112566. [Google Scholar] [CrossRef]

4. Miao X, Wang J, Bai H, Wang Y, Wu Z, Lin Q, et al. The numerical reactor system for China’s new-generation exascale supercomputer. Innovation. 2026;7(4):101268. doi:10.1016/j.xinn.2026.101268. [Google Scholar] [PubMed] [CrossRef]

5. Demaziere C, Sanchez-Espinoza VH, Zentner I. Advanced numerical simulation and modeling for reactor safety—contributions from the CORTEX, McSAFER, and METIS projects. EPJ Nuclear Sci Technol. 2022;8:29. doi:10.1051/epjn/2022026. [Google Scholar] [CrossRef]

6. Imron M. Development and verification of open reactor simulator ADPRES. Ann Nucl Energy. 2019;133(3):580–8. doi:10.1016/j.anucene.2019.06.049. [Google Scholar] [CrossRef]

7. Schneidesch C, Zhang J. Development and application of Tractebel’s multi-physics simulation capacity for advanced reactor safety evaluation. Nucl Eng Des. 2024;424:113253. doi:10.1016/j.nucengdes.2024.113253. [Google Scholar] [CrossRef]

8. Chan Y, Xiao S. Implementation and performance study of the lp-CMFD acceleration scheme for Monte Carlo method based k-eigenvalue neutron transport calculation in 1D geometry. Ann Nucl Energy. 2021;163(2):108562. doi:10.1016/j.anucene.2021.108562. [Google Scholar] [CrossRef]

9. Than YR, Xiao S. lp–CMFD acceleration schemes in multi-energy group 2D monte carlo transpor. Front Energy Res. 2022;10:1035797. doi:10.3389/fenrg.2022.1035797. [Google Scholar] [CrossRef]

10. Ong TKC, Xiao S, Peterson PF. An open-source thermo-hydraulic uniphase advection and convection solver for salt flows (TUAS). Int J Adv Nucl React Des Technol. 2024;6(4):281–301. doi:10.2139/ssrn.4998548. [Google Scholar] [CrossRef]

11. Finnemann H, Galati A. Benchmark problem book, report anl-7416 (suppl. 2). Argonne, IL, USA: Argonne National Laboratory; 1977. [Google Scholar]

12. Finnemann H, Galati A. NEACRP 3-D LWR core transient benchmark. Boulogne-Billancourt, France: OECD Nuclear Energy Agency; 1991. [Google Scholar]

13. Finnemann H, Bauer H, Galati A, Martinelli R. Results of LWR core transient benchmarks. Boulogne-Billancourt, France: OECD Nuclear Energy Agency; 1993. [Google Scholar]

14. Rubin A, Schoedel A, Avramova M, Utsuno H. OECD/NRC benchmark based on NUPEC PWR sub-channel and bundle tests (PSBT). Boulogne-Billancourt, France: OECD Nuclear Energy Agency; 2012. [Google Scholar]

15. Rubin A, Avramova M, Velazquez-Lozada A. International benchmark on pressurised water reactor sub-channel and bundle tests. Boulogne-Billancourt, France: OECD Nuclear Energy Agency; 2016. [Google Scholar]

16. Lawrence RD. Progress in nodal diffusion methods for the solution of the neutron diffusion and transport equations. Prog Nucl Energy. 1986;17(3):271–301. doi:10.1016/0149-1970(86)90034-x. [Google Scholar] [CrossRef]

17. Kim YI, Kim YJ, Kim SJ, Kim TK. A semi-analytic multigroup nodal method. Ann Nucl Energy. 1999;26(8):699–708. doi:10.1016/s0306-4549(98)00088-7. [Google Scholar] [CrossRef]

18. Fu XD, Cho NZ. Nonlinear analytic and semi-analytic nodal methods for multigroup neutron diffusion calculations. J Nucl Sci Technol. 2002;39(10):1015–25. doi:10.3327/jnst.39.1015. [Google Scholar] [CrossRef]

19. Tang C. Development and verification of an SP3 code using semi-analytic nodal method for pin-by-pin calculation. J Phys Sci Appl. 2017;7(2):108042. doi:10.17265/2159-5348/2017.02.002. [Google Scholar] [CrossRef]

20. R7-97(2012). Revised release on the IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam. Berlin, Germany: International Association for the Properties of Water and Steam; 2012. [Google Scholar]

21. Zou L, Zhao H, Zhang H, Lu Q. C++ implementation of IAWPS water/steam properties. Idaho Falls, ID, USA: Idaho National Laboratory; 2014. [Google Scholar]

22. Zou L, Zhao H, Zhang H. Numerical implementation, verification and validation of two-phase flow fourequation drift flux model with Jacobian-free Newton–Krylov method. Ann Nucl Energy. 2016;87(2):707–19. [Google Scholar]

23. NUREG/CR-5535. RELAP5/MOD3.3 code manual volume I. Rockville, MD, USA: Nuclear Regulatory Commission; 2001. [Google Scholar]

24. ML120060218. TRACE V5.0 theory manual. Norfolk, VA, USA: Commission USNR; 2010. [Google Scholar]

25. Hibiki T, Iishi W. One-dimensional drift-flux model and constitutive equations for relative motion between phases in various two-phase flow regimes. Int J Heat Mass Transf. 2003;46:4935–48. doi:10.2172/6871478. [Google Scholar] [CrossRef]

26. Hibiki T, Tsukamoto N. Drift-flux model for upward dispersed two-phase flows in vertical medium-to-large round tubes. Prog Nucl Energy. 2023;158:104611. doi:10.1016/j.pnucene.2023.104611. [Google Scholar] [CrossRef]

27. Rassame S, Hibiki T. Drift-flux model for dispersed adiabatic and boiling two-phase flows in rectangular channels. Int J Heat Mass Transf. 2024;224:125270. doi:10.1016/j.ijheatmasstransfer.2024.125270. [Google Scholar] [CrossRef]

28. Zhang H, Hibiki T, Xiao Y, Gu H. Two-group drift-flux model in tight lattice subchannel. Int Commun Heat Mass Transfer. 2024;159(C):108201. doi:10.1016/j.icheatmasstransfer.2024.108201. [Google Scholar] [CrossRef]

29. Yu M, Hibiki T. Two-group drift-flux model for dispersed gas-liquid flows in rod bundles. Int J Heat Mass Transf. 2024;222:125174. [Google Scholar]

30. Chexal B, Lellouche G. Full-range drift-flux correlation for vertical flows. Palo Alto, CA, USA: Electric Power Research Institute; 1985. [Google Scholar]

31. Vargaftik NB, Volkov BN, Voljak LD. Formula for water surface tension from international tables of the surface tension of water. J Phys Chem Ref. 1983;12(2):817–20. doi:10.1063/1.555688. [Google Scholar] [CrossRef]

32. Cheng NS. Formulas for friction factor in transitional regimes. J Hydraul Eng. 2008;134(9):1357–62. doi:10.1061/(asce)0733-9429(2008)134:9(1357). [Google Scholar] [CrossRef]

33. Todreas NE, Kazimi MS. Nuclear systems I—thermal hydraulic fundamentals. New York, NY, USA: Hemisphere Publishing Corporation; 1990. [Google Scholar]

34. Zhao X, Shirvan K, Wu Y, Kazimi MS. Critical power and void fraction prediction of tight bundle designs. Nucl Technol. 2016;196:553–67. [Google Scholar]

35. Eyert V. A comparative study on methods for convergence acceleration of iterative vector sequences. J Comput Phys. 1996;124(12):271–85. doi:10.1006/jcph.1996.0059. [Google Scholar] [CrossRef]

36. Anderson DG. Iterative procedures for nonlinear integral equations. J ACM. 1965;12(4):547–60. doi:10.1145/321296.321305. [Google Scholar] [CrossRef]

37. Facchini A, Lee J, Joo HG. Investigation of anderson acceleration in neutronics-thermal hydraulics coupled direct whole core calculation. Ann Nucl Energy. 2021;153:108042. [Google Scholar]

38. Knoll DA, Keyes DE. Jacobian-free newton-krylov methods: a survey of approaches and applications. J Comput Phys. 2004;193(2):357–97. [Google Scholar]

39. Islam A, Nushrat R, Rahim TA, Mollah AS. Modeling and validation of IAEA 3D PWR benchmark problem using COMSOL multiphysics code. Int J Integr Sci Technol. 2022;4:40–4. doi:10.1515/9781937585730-011. [Google Scholar] [PubMed] [CrossRef]

40. Joo HG, Barber D, Jiang G, Downar T. PARCS: a multi-dimensional two-group reactorkinetics code based on the nonlinear analytic nodal method. West Lafayette, IN, USA: Purdue University, School of Nuclear Engineering; 1998. [Google Scholar]

41. Tong LS. Prediction of departure from nucleate boiling for an axially non-uniform heat flux distribution. J Nucl Energy. 1967;21(3):241–8. doi:10.1016/s0022-3107(67)90054-8. [Google Scholar] [CrossRef]

42. Knight M. Derivation of a refined PANTHER solution to the NEACRP PWR rod ejection transients. In: Proceedings of the Joint International Conference on Mathematical Methods and Supercomputing in Nuclear Applications; 1997 Apr 19–23; Karlsruhe, Germany. [Google Scholar]

43. Guelfi A, Bestion D, Boucker M, Boudier P, Fillion P, Grandotto M, et al. NEPTUNE: a new software platform for advanced nuclear thermal hydraulics. Nucl Sci Eng. 2007;156(3):281–324. [Google Scholar]

44. Emonot P, Souyri A, Gandrille JL, Barré F. CATHARE-3: a new system code for thermal-hydraulics in the context of the NEPTUNE project. Nucl Eng Des. 2011;241(11):4476–81. [Google Scholar]

45. Toumi I, Bergeron A, Gallo D, Royer E, Caruge D. FLICA-4: a three-dimensional two-phase flow computer code with advanced numerical methods for nuclear applications. Nucl Eng Des. 2000;200(1–2):139–55. [Google Scholar]

46. Sulaiman IM, Mamat M, Omesa UA. Nonlinear systems–theoretical aspects and recent applications. London, UK: IntechOpen; 2020. [Google Scholar]

47. Yuan G, Hang X. Acceleration methods of nonlinear iteration for nonlinear parabolic equations. J Comput Math. 2006;24(3):412–24. [Google Scholar]

48. Vargun D. Acceleration methods for nonlinear solvers and application to fluid flow simulations. Clemson, SC, USA: Clemson University; 2023. [Google Scholar]

49. Radermacher A, Reese S. POD-based model reduction with empirical interpolation appliedto nonlinear elasticity. Int J Numer Meth Engng. 2016;107(6):477–95. doi:10.1002/nme.5177. [Google Scholar] [CrossRef]

50. Zhang Q, Ritzert S, Zhang J, Kehls J, Reese S, Brepols T. A unified multi-perspective quadratic manifold for mitigating the Kolmogorov barrier in multiphysics damage. J Mech Phys Solids. 2026;209(8):106499. doi:10.1016/j.jmps.2025.106499. [Google Scholar] [CrossRef]