Past nuclear accident occurrences raised strong concerns which led to research on nuclear safety. One of the major causes of nuclear accidents is the impeded circulation of core coolant, leading to decay heat removal cessation and rapid temperature rise. If uncontrolled, this results in critical heat flux, loss of coolant accidents, and core dryout. Detailed melted core relocation (i.e., nuclear fuel, graphite, and zircaloy) needs to be investigated through interface capture and multimaterial flow model coupling, which have not been done in previous studies. This work aims to investigate the impacts of temperature and core material composition on the flow dynamics during core relocation. In this study, mass fraction is discretized using a streamlined upwind Petrov–Galerkin method spatially and a modified Crank–Nicolson method temporally to accurately capture fluid interfaces using a high-order accurate flux-limiter. Two core material composition cases (individual material properties case and bulk material properties case) were considered to assess the impact of temperature and core materials composition on both flow dynamics and computational time. Temperature has a significant impact on core material transport and corium flow dynamics during core relocation. Bulk materials properties case has greater impact of temperature on its corium resulting in faster materials transport, but with higher computation time.

One strategy to address the global energy crisis is through nuclear power, despite the fact that its high heat output during operation makes it susceptible to catastrophic harm. In light of nuclear reactor accidents at Fukushima, Chernobyl, Three Mile Island, and other locations, risk and safety studies are vital to understand the origins and progression of such accidents.1,2 In nuclear water reactors, fuel rods contain nuclear material (i.e., enriched uranium fuel), which undergo fission reactions with release of heat. Cooling systems (with water as a coolant) extract this heat and maintain reactors at safe temperature. In pressurized water reactors (PWRs), if there is a loss of coolant due to a failure in the cooling system, the heat generated by the ongoing nuclear reactions may cause a significant rise in the temperature of the reactor core. This might result in core dryout, loss-of-coolant accidents (LOCA), and critical heat flux (CHF) within the coolant channel if not appropriately managed at an early stage.3 

Core materials, including fuel rods and structural components, are subjected to extreme thermal stress conditions. This results in changes of thermophysical properties, such as thermal conductivity, specific heat, and mechanical strength, leading to phase transitions, including melting. Molten fuel and cladding form the corium (a mixture of melted fuel, cladding, and reactor coolant), which, being highly radioactive, can potentially breach the containment vessel. Core relocation occurs as a result of the interaction between corium and structural materials, which compromises the integrity of the reactor containment.

Fuel rods overheat during CHF events, triggering LOCA, which is associated with a sharp reduction in local heat flux as a result of nucleation boiling. Common features of LOCA include rapid system pressure drop, abrupt local heat flux reduction, core dryout, cladding and fuel rod ballooning, core relocation, high temperature steam oxidation, hydrogen release, and core degradation.3–5 LOCA leads to reduced reactor's cooling capacity with large heat accumulation, resulting in local hydrodynamic instabilities.6,7 In assessing CHF, post-CHF, LOCA events, and core relocation, multiphase interface capturing modeling of the reactor core may help predict flow dynamics in the reactor core. Modeling the behavior of the material interface during severe accidents is essential for fundamental understanding of material interactions and system flow dynamics.8,9

In order to present an adequate analysis of the multiphase flow dynamics during core relocation events, it is crucial to track the interfaces between different materials (i.e., melted nuclear fuel and cladding). Interface modeling primarily employs two methods: Eulerian and Lagrangian tracking techniques. The element (mesh cell) carrying discretized quantities is transported over the continuous media in Lagrangian schemes, whereas the grid is fixed in the Eulerian approach, while the fluid moves between computational cells.10 Numerous fluid dynamics and transport problems have been solved using interface tracking and capturing techniques.11–14 Some of the main types of interface capture techniques are volume of fluid (VOF), level set (LS), and phase field methods. Among these, VOF is the most efficient and reliable approach due to its mass conservation ability across the interfaces.14 Using an implicit technique to adjust mesh refinement during calculations with accurate mass conservation, VOF methods create moving contact lines along no-slip borders.15 With every successive iteration in VOF, the interface between two fluids is reconstructed using spatial information.

Advanced volume fraction expansion functions schemes, such as piece-wise linear interface calculation (PLIC), coupled level set and volume of fluid (CLS-VOF) method, and simplified VOF (SVOF), were developed to overcome interface curvature challenges encountered when using the classic VOF scheme.10,16,17 An interface capture approach based on a switching strategy had been developed by Darwish and Moukalled18 with the Switching Technique for Advection and Capturing of Surfaces (STACS) scheme. Compressive advection interface capturing method (CAICM), a special VOF variant, has been applied to interface capturing of collapsing water columns by Pavlidis et al.19 In order to develop bounded and accurate interface capture solutions, this computational methodology was developed using a high-order control volume finite element method (CVFEM) to obtain fluxes subject to flux-limiting (based on a normalized variable diagram) on control volume (CV) boundaries.

This research aims to investigate the impacts of temperature and core material composition on the flow dynamics during core relocation. The simulations were conducted using CAICM coupled with high-order accurate schemes discretized with the streamline upwind Petrov–Galerkin (SUPG) method.

The remainder of this paper is structured as follows: Secs. II and III describe the model formulation and validation, respectively, while Secs. IV and V discuss core relocation and concluding remarks, respectively.

The model formulation in this study utilizes higher-order CVFEM and CAICM with flux limiters using a normalized variable diagram approach.

The formulation is based on a family of unstructured triangles (2D) and tetrahedral (3D) element types, PnDG-Pm and PnDG-PmDG.20 Each finite element is split up into a number of CVs, which can be quadrilaterals (2D) or hexahedra (3D) as shown in Fig. 1. For these element types, pressure is represented by a mth-order polynomial in either continuous or discontinuous space, whereas velocity is represented by a nth-order polynomial in discontinuous space. In this study, P1DG-P2 (Fig. 1 linear discontinuous velocity function and quadratic continuous pressure function) elements were used for the model discretization. In this formulation, scalars (such as temperature, density, volume fraction, etc.) are discretized in CV space, whereas pressure and velocity are discretized in FE space.

FIG. 1.
Sketch of P1DGP2 element type. The tick pink line triangles are finite elements, quadrilateral dotted lines are control volumes, and the black circles represent quadratic continuous pressure nodes, while the dimmed-red circles represent linear discontinuous velocity nodes.
View largeDownload slide

Sketch of P1DGP2 element type. The tick pink line triangles are finite elements, quadrilateral dotted lines are control volumes, and the black circles represent quadratic continuous pressure nodes, while the dimmed-red circles represent linear discontinuous velocity nodes.

FIG. 1.
Sketch of P1DGP2 element type. The tick pink line triangles are finite elements, quadrilateral dotted lines are control volumes, and the black circles represent quadratic continuous pressure nodes, while the dimmed-red circles represent linear discontinuous velocity nodes.
View largeDownload slide

Sketch of P1DGP2 element type. The tick pink line triangles are finite elements, quadrilateral dotted lines are control volumes, and the black circles represent quadratic continuous pressure nodes, while the dimmed-red circles represent linear discontinuous velocity nodes.

Close modal
The model formulation is based on a multimaterial approach, in which the system is represented as a single phase made up of numerous immiscible materials (fluids) separated by distinct interfaces. This formulation models materials with discrete transported scalars (e.g., temperature and mass fraction) but with single velocity and pressure. Conservation equations were solved for the multimaterial components, but only the transport equation is shown here. Multimaterial formulations for continuity and momentum equations are fully described in Refs. 19, 21, and 22. In the absence of a source term, the advection-diffusion transport is expressed as
(1)
where α is the transported scalar quantity [i.e., mass fraction and enthalpy (CpdT)] and κ is the diffusion term. u → is the flow velocity vector, ρ is the material density, Cp is the heat capacity, T is the temperature, and t is the time. Bulk thermophysical properties of the fluids can be expressed as
(2)
where αj, ρj, κj, and C p , j represent the mass fraction, density, diffusion term, and heat capacity, respectively, of the material j. j = 1 , 2 , … , N⁠, where N represents the total number of materials.

Discontinuities of flow fields across the fluid interfaces, complex topologies, and scale separation around CV surfaces are addressed using the Petrov–Galerkin spatial discretization approach. A modified Crank–Nicolson scheme with a time-stepping parameter τ is used to discretize Eq. (1) temporally and compute the advected quantity within the CV.

With the assumption of no source term in Eq. (1), the transport equation is solved with the introduction of CV shape functions Mi. With Green's theorem, Eq. (1) can be simplified using the Petrov–Galerkin spatial discretization approach and integration by parts to give Eq. (3),
(3)
The stiffness matrix from the third term of the integral in Eq. (3) is symmetrical since α and Mi are exchangeable without any effects on the integral, whereas that of the second term is nonsymmetrical. Equation (3) is further simplified to
(4)
where ∂ V i is the outward normal surface area of the control volume, and i = 1 , 2 , … , I⁠, where I is the number of control volumes. n is the outward pointing unit vector to the surface (⁠ ∂ V i⁠) of the control volume. The difference between the advected and diffusive quantity within the FE space in the last term of Eq. (4) gives the transported quantity, α *⁠, across the interface as shown in Fig. 2, represented in 1D. Introducing a time-stepping parameter τ of the value 1/2, Eq. (4) can be applied over each control volume, i, as
(5)
To compute the fluxes across the interface of node associated with CV i, extra nodes within CV, i − 1 2 and i + 1 2⁠, were added between i and i − 1, and i and i + 1, respectively, as shown in Fig. 2. For a particular CV i with the edge length of Δ x (⁠ x i − x i + 1⁠), the advected quantity α across the interface of node between i − 1 2 and i + 1 2 at time step k is set as α i + 1 2 k and α i − 1 2 k⁠, respectively, as shown in Fig. 2. The discretized form of Eq. (1) becomes
(6)
where ( n · u i + 1 2 k + 1 α j , i + 1 2 k + 1 − n · ϑ j , i + 1 2 k + 1 ∇ α k + 1 | j , i + 1 2 ) =  q j , i + 1 2 k + 1 and ( n · ϑ j , i + 1 / 2 k ∇ α k | i + 1 2 − n · u i + 1 2 k α j , i + 1 2 k ) = − q j , i + 1 2 k are flux-limiting values on the face i + 1 2 at discrete times k + 1 and k, respectively. ∇ α k + 1 | j , i + 1 2 and ∇ α k | j , i + 1 2 are the gradient of the transported material at face i + 1 2 for k + 1 and k time steps, respectively, for material j. ϑ = κ ρ⁠. The value of τ for face i − 1 2 and face i + 1 2 at time step from k to k + 1 is given by τ j , i − 1 2 k + 1 2 and τ j , i + 1 2 k + 1 2⁠, respectively. Equation (6) results in the transported quantities evaluated on the CV i at the time step k + 1 relative to surrounding faces i + 1 2 and i − 1 2 as a function of their flux-limiting values. n can be expressed as the vector space normalization in terms of the transported scalar gradient ∇ α ( n = ∇ α | ∇ α | )⁠. Transported quantity is then given as
(7)
The flux-limiting scheme was applied to the scalar's and flux transport across the interface of CV i connecting upwind and downwind CVs. Further details of the flux-limiter application with normalized variable diagram (NVD) approach are described in Refs. 19 and 23.
FIG. 2.
1D finite element diagram with downwind (wd), upwind (wd), and centered (g) CVs, with higher-order FE field reconstruction face value ( α *).
View largeDownload slide

1D finite element diagram with downwind (wd), upwind (wd), and centered (g) CVs, with higher-order FE field reconstruction face value (⁠ α *⁠).

FIG. 2.
1D finite element diagram with downwind (wd), upwind (wd), and centered (g) CVs, with higher-order FE field reconstruction face value ( α *).
View largeDownload slide

1D finite element diagram with downwind (wd), upwind (wd), and centered (g) CVs, with higher-order FE field reconstruction face value (⁠ α *⁠).

Close modal

For the validation of this model, the classical collapsing water column (CWC) was used based on experimental and simulation work by Cruchaga et al.24 and Yeoh and Barber,25 respectively. Cruchaga et al.24 designed an experiment with a water column height H to length L ratio of 2 (i.e., Ar = 2, A r = H L⁠). The former25 used the same water column height H to length L ratio (as used in this study, Ar = 2) for numerical interface capturing investigations using FVM-based ANSYS, Inc., CFX model. Water column of height H was held with a barrier at one end (left side) of a tank. The time for the water column collapsing at different vertical and horizontal positions of the tank was recorded as the barrier was removed, which was simulated numerically. The laboratory experiment was performed such that the time evolution of the collapsing column was captured by a video camera and stopwatch installed on the tank. A qualitative validation is presented to compare the present study against the lab-scale and simulation experiments. The coupled CVFEM and CAICM model formulation was initially validated by Pavlidis et al.19 (simple advection and water collapsing column with multiphase flow model and anisotropic unstructured mesh), Xie et al.21 (falling liquid film, rising bubble, milkcrown phenomena with single phase, and multi-component formulation), and Bregu et al.5 (water collapsing column, single/multiphase, and multi-component formulations).

The collapsing of the water column was under the influence of gravitational force. The physics behind collapsing water column and core relocation is the same. At 0.2 s, the tank floor has almost been covered by the collapsed water in the three cases as seen in Fig. 3. The water column collapsed under the influence of gravitational force. From Fig. 3, it can be observed that the current study model is in good agreement with both studies.

FIG. 3.
Model validation using cwc simulation experiments at 0.2 s time for the present study (left), Yeoh and Barber25 (center), and Cruchaga et al.24 (right).
View largeDownload slide

Model validation using cwc simulation experiments at 0.2 s time for the present study (left), Yeoh and Barber25 (center), and Cruchaga et al.24 (right).

FIG. 3.
Model validation using cwc simulation experiments at 0.2 s time for the present study (left), Yeoh and Barber25 (center), and Cruchaga et al.24 (right).
View largeDownload slide

Model validation using cwc simulation experiments at 0.2 s time for the present study (left), Yeoh and Barber25 (center), and Cruchaga et al.24 (right).

Close modal

When the reactor core is subjected to extreme thermal conditions, as explained in Sec. I, core materials can undergo melting. This melting occurs as a result of a change in the thermophysical properties of the materials. Investigating the behavior of nuclear reactor cores under extreme conditions, such as severe accidents or core relocation, is vital for the reactor safety improvement. Assessing the possibility and mechanisms of melted core materials relocating within the reactor vessel or containment structure during such events could help in maintaining reactor integrity. Exploring some factors such as heat flow, melted materials, and structural components interactions, and potential consequences of such relocation on containment integrity are significant to nuclear safety improvement. Using models and simulations to predict the behavior of melted materials during core relocation can ultimately inform safety measures, emergency response strategies, and reactor design improvements to mitigate severe accidents consequences.

In this study, the impact of corium composition on the flow dynamics during core relocation is numerically investigated. Simulations involving 2 case scenarios were conducted as shown in Table I: The containment (including core material) in scenario 1 is comprised of UO2 and air, while that of scenario 2 comprised of the combination of UO2, zircaloy, graphite, water (of 10%, 10%, 15%, 65% composition, respectively), and air. Mesh resolution, mean flow velocity, and temperature of the two scenarios are investigated. This is necessary to understand the response of different corium composition and impacts of temperature on the flow dynamics during core relocation. Figure 4 shows the initial unstructured mesh resolution and material configuration (at t = 0 s) for a grid-independent study. The domain is a 2D geometry of L = 1   m (see Fig. 4) with 12 050 triangular elements. Melted core materials are in the red-colored region, while the rest of the domain is filled with air (in blue) at the onset of the core relocation. The corium collapse may be presented as symmetrical (as the corium is assumed with homogeneous thermophysical and rheological properties); however, corium motion after it hits the bottom is not symmetrical. Here, we assume symmetry of the geometry, and the next research stage will be a full 3D simulation (future work).

TABLE I.

Descriptions of the different scenarios of core relocation simulations.

Descriptions Considered scenarios
Scenario 1 Scenario 2
Core materials composition  UO2 + air  (water, UO2, graphite and zircaloy) + air 
Core materials properties  Individual material property  Bulk materials properties 
Percentage composition  100% UO2 + air  [water (65%), UO2 (10%), graphite (15%), zircaloy (10%)] + air 
Descriptions Considered scenarios
Scenario 1 Scenario 2
Core materials composition  UO2 + air  (water, UO2, graphite and zircaloy) + air 
Core materials properties  Individual material property  Bulk materials properties 
Percentage composition  100% UO2 + air  [water (65%), UO2 (10%), graphite (15%), zircaloy (10%)] + air 
View Large
FIG. 4.
Volume fraction representation of melted core with unstructured adaptive mesh. There are (initially) 12 050 triangular P1DG-P2 elements.
View largeDownload slide

Volume fraction representation of melted core with unstructured adaptive mesh. There are (initially) 12 050 triangular P1DG-P2 elements.

FIG. 4.
Volume fraction representation of melted core with unstructured adaptive mesh. There are (initially) 12 050 triangular P1DG-P2 elements.
View largeDownload slide

Volume fraction representation of melted core with unstructured adaptive mesh. There are (initially) 12 050 triangular P1DG-P2 elements.

Close modal
In the simulations presented here, core materials (UO2, graphite, water, and zircaloy) are assumed incompressible, whereas air is assumed to behave according to the following equation of state:
(8)
where ρ ° is the reference density and a and b are elasticity bulk modulus and thermal expansion coefficient given as 7.06 × 10 − 6   Pa − 1 and 9.4 × 10 − 4   K − 1⁠, respectively. P and T are pressure and temperature of the system, respectively, while T ° and P ° are the reference temperature (323 K) and pressure (1.0  × 10 5 Pa), respectively.

Temperature Dirichlet boundary condition of 1000 K is imposed on the walls of the reactor. Core materials were assumed to be fully melted (initially at a temperature of 2700 K) with the air temperature of 1200 K and collapsed at the same time with no nuclear recriticality. The reactor coolant (water) was also assumed to exist in liquid form throughout the relocation event.

An anisotropic adaptive mesh method as described by Piggott et al.26 was used to optimally control mesh resolution based on the flow dynamics. Adaptive time steps with Courant–Friedrichs–Lewy (CFL) number of 5.0 were applied to simulations to capture sharp interfaces between core materials with good convergence. The method is able to dynamically focus mesh resolution at the interface between air and corium, whereas in the remaining domain, a coarse mesh can be used, reducing computational cost. A detailed description of the mesh adaptivity method used in this work can be found.27–29 In the simulations shown in this section, the mesh is adapted based on the UO2 mass fraction with an interpolation error bound of 0.0005, a gradation parameter of 1.5, an aspect ratio bound of 3, and the mesh was adapted at every three time steps. Simulations were conducted using Intel Core i5 with 8 cores and 16 GB RAM. Simulations of scenario 1 required 84 h of wall clock time to complete. For scenario 2, simulations required 94.5 h of wall clock time. This shows the higher computation cost of the bulk material simulation. Simulations for both case scenarios were initialized with mesh resolution of 12 050 triangular P1DGP2 element types.

Figure 5 shows time evolution of the number of CVs and elements for both cases (scenarios 1 and 2). Due to mesh resolution, the number of elements and CVs increase gradually at the early stage of the simulation. Mesh resolution follows similar pattern for both cases, though with larger number of elements used in scenario 2. There was a spike in number elements and CV evolution as the corium materials approach and hit the lower plenum (Fig. 6). This shows slightly faster material relocation of scenario 1 than scenario 2 for the corresponding relocation time.

FIG. 5.
Number of control volumes (CVs) and elements evolution. A and B represent number of elements and CVs for scenario 1, while C and D represent number of elements and CVs for scenario 2.
View largeDownload slide

Number of control volumes (CVs) and elements evolution. A and B represent number of elements and CVs for scenario 1, while C and D represent number of elements and CVs for scenario 2.

FIG. 5.
Number of control volumes (CVs) and elements evolution. A and B represent number of elements and CVs for scenario 1, while C and D represent number of elements and CVs for scenario 2.
View largeDownload slide

Number of control volumes (CVs) and elements evolution. A and B represent number of elements and CVs for scenario 1, while C and D represent number of elements and CVs for scenario 2.

Close modal
FIG. 6.
Scenario 1 adaptive mesh at t = 0.4 s (left) and 0.99 s, and mass fraction distribution at t = 1.40 s (right).
View largeDownload slide

Scenario 1 adaptive mesh at t = 0.4 s (left) and 0.99 s, and mass fraction distribution at t = 1.40 s (right).

FIG. 6.
Scenario 1 adaptive mesh at t = 0.4 s (left) and 0.99 s, and mass fraction distribution at t = 1.40 s (right).
View largeDownload slide

Scenario 1 adaptive mesh at t = 0.4 s (left) and 0.99 s, and mass fraction distribution at t = 1.40 s (right).

Close modal

Figure 7 shows mean flow velocities for both case scenarios. There was an increase in the velocity at 0.70 s as the corium hits the bottom of the plenum due to the splashing of the materials as shown in Fig. 6. A larger spike in the mean velocity was observed for the simulation conducted for scenario 2. Simulations demonstrated that corium in scenario 1 reaches the lower plenum slightly faster than in scenario 2 with max mean velocity for both scenarios at 0.70 s with a small (negligible) time-lag between them. This is expected as density (or weight) does not affect the time of free fall. Air resistance has limited/small impact on the speed of the fall.

FIG. 7.
Mean flow velocity for scenarios 1 (individual) and 2 (bulk).
View largeDownload slide

Mean flow velocity for scenarios 1 (individual) and 2 (bulk).

FIG. 7.
Mean flow velocity for scenarios 1 (individual) and 2 (bulk).
View largeDownload slide

Mean flow velocity for scenarios 1 (individual) and 2 (bulk).

Close modal

The mean temperature over simulation time for both case scenarios is shown in Fig. 8. Material 1 (UO2 for scenario 1 and [water (65% w/w), UO2 (10% w/w), graphite (15% w/w), zircaloy (10% w/w)] for scenario 2) was assumed to have melted and mix with air while collapsing. A mean temperature fluctuation was observed during the relocation. Heat exchange between air and corium leads to mean temperature ranging from 1550 to 2265 K. There was a higher temperature recorded for the bulk material case. This is attributed to the larger effective thermal conductivity in scenario 2 than in scenario 1 [ κ eff scenario   1 = 6.80 W   ( m   K ) − 1     and     κ eff scenario   2 = 23.424 W   ( m   K ) − 1⁠, Table II].

FIG. 8.
Corium and air mean temperature for both cases.
View largeDownload slide

Corium and air mean temperature for both cases.

FIG. 8.
Corium and air mean temperature for both cases.
View largeDownload slide

Corium and air mean temperature for both cases.

Close modal
TABLE II.

Thermophysical properties of air and corium materials at 2123 K and 100 bar.30,31

Properties Core materials percentage compositions
Core materials Water UO2 Zircaloy Graphite Air
Materials % composition (w/w)  65%  10%  10%  15%   
Density (⁠ kgm − 3⁠)  958.78  19 100  6560  1710  1.2945 (reference density) 
Specific heat capacity [ J ( kg   K ) − 1]  4184  280  285  720  1006 
Thermal conductivity [ W ( m   K ) − 1]  0.591  6.80  13.6  140  0.025 
Viscosity [ kg ( m   s ) − 1]  2.822 × 10 − 4  6.5  10.3  10  2.181 × 10 − 5 
Properties Core materials percentage compositions
Core materials Water UO2 Zircaloy Graphite Air
Materials % composition (w/w)  65%  10%  10%  15%   
Density (⁠ kgm − 3⁠)  958.78  19 100  6560  1710  1.2945 (reference density) 
Specific heat capacity [ J ( kg   K ) − 1]  4184  280  285  720  1006 
Thermal conductivity [ W ( m   K ) − 1]  0.591  6.80  13.6  140  0.025 
Viscosity [ kg ( m   s ) − 1]  2.822 × 10 − 4  6.5  10.3  10  2.181 × 10 − 5 
View Large

In this study, impacts of temperature and core material composition on the flow dynamics during post-CHF events in PWRs are numerically investigated. Two corium material compositions (UO2 and air for scenario 1 and [water (65% w/w), UO2 (10% w/w), graphite (15% w/w), zircaloy (10% w/w)] and air for scenario 2) were considered to compare the impacts of temperature and their computation costs. The interface between corium materials is modeled through an integrated high-order accurate control volume finite element method (CVFEM) and the compressive advection interface capturing method (CAICM). To control numerical instabilities, the model was discretized using SUPG and the flux-limiting schemes with a normalized variable diagram approach. Mesh adaptivity algorithm was applied to focus resolution at the corium–air interfaces. Using EoS for the gas phase, the effects of air density variation on the corium mass fraction transport were determined. From the results of this study, the following conclusions are made:

  • Air compressibility has a greater impact on the transport of scenario 2 corium than that of scenario 1 during core relocation events.

  • Splashing of the corium after hitting the bottom of the plenum resulted in mean flow velocity increase and temperature fluctuation.

  • Simulated mean corium temperature and velocity have similar orders of magnitude for both scenarios.

In summary, based on the results presented in Sec. IV, fluid and thermal energy flow dynamics (mean flow velocity and temperature) are similar for scenarios 1 (i.e., lumping all materials into UO2) and 2 (explicitly defining properties of all corium materials). However, computational cost for simulating scenario 1 is significantly smaller than for scenario 2 (84 h and 94.5 h of wall-time processed in an Intel Core i5 with 8 cores and 16 GB RAM). Multiphase interface capturing requires high mesh resolution and computational cost. Mesh adaptivity is used to focus grid resolution at the relevant physics, in this case the interface between all materials of the corium, to minimize computational cost without compromising numerical accuracy. As a good CFD practice, we start simulations with mesh-grids that is deemed as grid-independent. Pain et al.,32 Hiester et al.,29 Adam et al.,28 and Xu et al.27 described how finite elements' shapes and topology change according to flow criteria to ensure numerical stability and convergence with controlled impact on numerical accuracy. Future work will extend the model to simulate core relocation in 3D.

The first author acknowledges the Petroleum Trust Development Fund (PTDF), Nigeria, for providing the funding for this research with Grant No. PTDF/ED/OSS/PHD/SAA/1801/20-20PHD109.

The authors have no conflicts to disclose.

Stephen A. Ajah: Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – original draft (equal). Lateef Temitope Akanji: Supervision (equal); Writing – review & editing (equal). Jefferson Gomes: Conceptualization (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
H.
Madokoro
 and
I.
Sato
, “
Estimation of the core degradation and relocation at the Fukushima Daiichi Nuclear Power Station Unit 2 based on RELAP/SCDAPSIM analysis
,”
Nucl. Eng. Des.
376
,
111123
(
2021
).
https://doi.org/10.1016/j.nucengdes.2021.111123
Google Scholar
Crossref
Search ADS
 
2.
S. J.
Hong
,
D. K.
Ko
,
G. H.
Lee
,
Y. A.
Kim
,
K. S.
Park
,
H. W.
Kim
,
C. G.
Kim
, and
G. H.
Jang
, “
An assessment of radiation exposure of workers during decommissioning of reactor vessel at Kori unit 1
,”
Nucl. Eng. Des.
392
,
111771
(
2022
).
https://doi.org/10.1016/j.nucengdes.2022.111771
Google Scholar
Crossref
Search ADS
 
3.
N.
Zhao
,
Y.
Chen
,
W.
Ma
,
S.
Bechta
, and
P.
Isaksson
, “
A nodal sensitivity study of MELCOR simulation for severe accidents in a pressurized water reactor
,”
Ann. Nucl. Energy
160
,
108373
(
2021
).
https://doi.org/10.1016/j.anucene.2021.108373
Google Scholar
Crossref
Search ADS
 
4.
S. H.
Gyepi-Garbrah
 and
T.
Nitheanandan
, “
Severe accident prevention and mitigation in pressurized heavy water reactors
,” in
Pressurized Heavy Water Reactors
(
Elsevier
,
2022
), pp.
477
–
508
.
https://doi.org/10.1016/B978-0-12-822054-2.00017-9
Google Scholar
Crossref
Search ADS
 
5.
E.
Bregu
,
S. A.
Ajah
, and
J.
Gomes
, “
Nuclear meltdown relocation and core catcher analysis
,”
VETOR
33
(
2
),
2
–
10
(
2023
).
https://doi.org/10.14295/vetor.v33i2.15158
Google Scholar
Crossref
Search ADS
 
6.
T. B.
Nguyen
,
D.
Liu
,
M. I.
Kayes
,
B.
Wang
,
N.
Rashin
,
P. W.
Leu
, and
T.
Tran
, “
Critical heat flux enhancement in pool boiling through increased rewetting on nanopillar array surfaces
,”
Sci. Rep.
8
,
4815
(
2018
).
https://doi.org/10.1038/s41598-018-22693-z
Google Scholar
Crossref
Search ADS
PubMed
 
7.
P.
Fernández-Arias
,
D.
Vergara
, and
J. A.
Orosa
, “
A global review of PWR nuclear power plants
,”
Appl. Sci.
10
,
4434
(
2020
).
https://doi.org/10.3390/app10134434
Google Scholar
Crossref
Search ADS
 
8.
J.
Wang
,
C.
Wang
,
K.
Shi
, and
G.
Su
, “
Review on core degradation and material migration research in light-water reactors
,”
Front. Energy Res.
6
,
1
–
10
(
2018
).
https://doi.org/10.3389/fenrg.2018.00003
Google Scholar
Crossref
Search ADS
 
9.
D.
Zhan
,
X.
Zhao
,
S.
Xia
,
P.
Chen
, and
H.
Chen
, “
Numerical simulation and validation for early core degradation phase under severe accidents
,”
Sci. Technol. Nucl. Install.
2020
,
6798738
.
https://doi.org/10.1155/2020/6798738
Crossref
Search ADS
 
10.
B. M.
Ningegowda
 and
B.
Premachandran
, “
A coupled level set and volume of fluid method with multi-directional advection algorithms for two-phase flows with and without phase change
,”
Int. J. Heat Mass Transfer
79
,
532
–
550
(
2014
).
https://doi.org/10.1016/j.ijheatmasstransfer.2014.08.039
Google Scholar
Crossref
Search ADS
 
11.
S.
McKee
,
M. F.
Tomé
,
V. G.
Ferreira
,
J. A.
Cuminato
,
A.
Castelo
,
F. S.
Sousa
, and
N.
Mangiavacchi
, “
The MAC method
,”
Comput. Fluids
37
,
907
–
930
(
2008
).
https://doi.org/10.1016/j.compfluid.2007.10.006
Google Scholar
Crossref
Search ADS
 
12.
S.
Noro
,
H.
Koibuchi
,
S.
Hongo
,
S.
Nagahiro
,
H.
Ikai
,
M.
Nakayama
,
T.
Uchimoto
, and
J.-P.
Rieu
, “
Langevin Navier-Stokes simulation of protoplasmic streaming by 2D MAC method
,” arXiv:2112.10901 (
2021
).
Google Scholar
 
13.
S.
Shin
,
J.
Chergui
,
D.
Juric
,
L.
Kahouadji
,
O. K.
Matar
,
R. V.
Craster
,
S.
Shin
,
J.
Chergui
,
D.
Juric
,
L.
Kahouadji
, and
O. K.
Matar
, “
An interface-tracking technique for multiphase flow with soluble surfacant
,”
J. Comput. Phys.
359
,
409
–
435
(
2018
).
https://doi.org/10.1016/j.jcp.2018.01.010
Google Scholar
Crossref
Search ADS
 
14.
Z.
Huang
,
G.
Lin
, and
A. M.
Ardekani
, “
A consistent and conservative phase-field model for thermo-gas-liquid-solid flows including liquid-solid phase change
,”
J. Comput. Phys.
449
,
110795
(
2022
).
https://doi.org/10.1016/j.jcp.2021.110795
Google Scholar
Crossref
Search ADS
 
15.
D. A.
Hoang
,
V.
van Steijn
,
L. M.
Portela
,
M. T.
Kreutzer
, and
C. R.
Kleijn
, “
Benchmark numerical simulations of segmented two-phase flows in microchannels using the volume of fluid method
,”
Comput. Fluids
86
,
28
–
36
(
2013
).
https://doi.org/10.1016/j.compfluid.2013.06.024
Google Scholar
Crossref
Search ADS
 
16.
I.
Chakraborty
,
G.
Biswas
, and
P. S.
Ghoshdastidar
, “
A coupled level-set and volume-of-fluid method for the buoyant rise of gas bubbles in liquids
,”
Int. J. Heat Mass Transfer
58
,
240
–
259
(
2013
).
https://doi.org/10.1016/j.ijheatmasstransfer.2012.11.027
Google Scholar
Crossref
Search ADS
 
17.
F.
Xiao
,
S.
Ii
, and
C.
Chen
, “
Revisit to the THINC scheme: A simple algebraic VOF algorithm
,”
J. Comput. Phys.
230
,
7086
–
7092
(
2011
).
https://doi.org/10.1016/j.jcp.2011.06.012
Google Scholar
Crossref
Search ADS
 
18.
M.
Darwish
 and
F.
Moukalled
, “
Convective schemes for capturing interfaces of free-surface flows on unstructured grids
,”
Numer. Heat Transfer, Part B
49
,
19
–
42
(
2006
).
https://doi.org/10.1080/10407790500272137
Google Scholar
Crossref
Search ADS
 
19.
D.
Pavlidis
,
J. L. M. A.
Gomes
,
Z.
Xie
,
J. R.
Percival
,
C. C.
Pain
, and
O. K.
Matar
, “
Compressive advection and multi-component methods for interface-capturing
,”
Numer. Methods Fluids
80
,
256
–
282
(
2016
).
https://doi.org/10.1002/fld.4078
Google Scholar
Crossref
Search ADS
 
20.
P.
Salinas
,
J. R.
Percival
,
D.
Pavlidis
,
Z.
Xie
,
J.
Gomes
,
C. C.
Pain
, and
M. D.
Jackson
, “
A discontinuous overlapping control volume finite element method for multi-phase porous media flow using dynamic unstructured mesh optimization
,” in
SPE Reservoir Simulation Symposium
(
Society of Petroleum Engineers
,
2015
), pp.
1511
–
1528
.
Google Scholar
Crossref
Search ADS
 
21.
Z.
Xie
,
D.
Pavlidis
,
P.
Salinas
,
J. R.
Percival
,
C. C.
Pain
, and
O. K.
Matar
, “
A balanced-force control volume finite element method for interfacial flows with surface tension using adaptive anisotropic unstructured meshes
,”
Comput. Fluids
138
,
38
–
50
(
2016
).
https://doi.org/10.1016/j.compfluid.2016.08.005
Google Scholar
Crossref
Search ADS
 
22.
H.
Lee
 and
S. H.
Rhee
, “
A dynamic interface compression method for VOF simulations
,”
J. Mech. Sci. Technol.
29
,
1849
–
1857
(
2015
).
https://doi.org/10.1007/s12206-015-0405-6
Google Scholar
Crossref
Search ADS
 
23.
M.
Darwish
 and
F.
Moukalled
, “
The χ schemes: A new consistent high-resolution formulation based on the normalized variable methodology
,”
Comput. Methods Appl. Mech. Eng.
192
(13–14),
1711
(
2003
).
https://doi.org/10.1016/S0045-7825(03)00199-3
Google Scholar
Crossref
Search ADS
 
24.
M. A.
Cruchaga
,
D. J.
Celentano
, and
T. E.
Tezduyar
, “
Collapse of a liquid column: Numerical simulation and experimental validation
,”
Comput. Mech.
39
,
453
–
476
(
2007
).
https://doi.org/10.1007/s00466-006-0043-z
Google Scholar
Crossref
Search ADS
 
25.
G. H.
Yeoh
 and
T.
Barber
, “
Assessment of interface capturing methods in Computational Fluid Dynamics (CFD) codes—A case study
,”
J. Comput. Multiphase Flows
1
,
201
–
215
(
2009
).
https://doi.org/10.1260/175748209789563946
Google Scholar
Crossref
Search ADS
 
26.
M. D.
Piggott
,
P. E.
Farrell
,
C. R.
Wilson
,
G. J.
Gorman
, and
C. C.
Pain
, “
Anisotropic mesh adaptivity for multi-scale ocean modelling
,”
Philos. Trans. R Soc. A
367
,
4591
–
4611
(
2009
).
https://doi.org/10.1098/rsta.2009.0155
Google Scholar
Crossref
Search ADS
 
27.
X.
Xu
,
C. C.
Pain
,
A. J.
Goddard
, and
C. R.
De Oliveira
, “
An automatic adaptive meshing technique for Delaunay triangulations
,”
Comput. Methods Appl. Mech. Eng.
161
,
297
–
303
(
1998
).
https://doi.org/10.1016/S0045-7825(97)00322-8
Google Scholar
Crossref
Search ADS
 
28.
A.
Adam
,
D.
Pavlidis
,
J. R.
Percival
,
P.
Salinas
,
Z.
Xie
,
F.
Fang
,
C. C.
Pain
,
A. H.
Muggeridge
, and
M. D.
Jackson
, “
Higher-order conservative interpolation between control-volume meshes: Application to advection and multiphase flow problems with dynamic mesh adaptivity
,”
J. Comput. Phys.
321
,
512
–
531
(
2016
).
https://doi.org/10.1016/j.jcp.2016.05.058
Google Scholar
Crossref
Search ADS
 
29.
H. R.
Hiester
,
M. D.
Piggott
,
P. E.
Farrell
, and
P. A.
Allison
, “
Assessment of spurious mixing in adaptive mesh simulations of the two-dimensional lock-exchange
,”
Ocean Modell.
73
,
30
–
44
(
2014
).
https://doi.org/10.1016/j.ocemod.2013.10.003
Google Scholar
Crossref
Search ADS
 
30.
L.
Ma
,
J.
Wang
,
A. M.
Marconnet
,
A. C.
Barbati
,
G. H.
McKinley
,
W.
Liu
, and
G.
Chen
, “
Viscosity and thermal conductivity of stable graphite suspensions near percolation
,”
Nano Lett.
15
,
127
–
133
(
2015
).
https://doi.org/10.1021/nl503181w
Google Scholar
Crossref
Search ADS
PubMed
 
31.
T.
Kondo
,
H.
Muta
,
K.
Kurosaki
,
F.
Kargl
,
A.
Yamaji
,
M.
Furuya
, and
Y.
Ohishi
, “
Density and viscosity of liquid ZrO2 measured by aerodynamic levitation technique
,”
Heliyon
5
,
e02049
(
2019
).
https://doi.org/10.1016/j.heliyon.2019.e02049
Google Scholar
Crossref
Search ADS
PubMed
 
32.
C. C.
Pain
,
A. P.
Umpleby
,
C. R.
de Oliveira
, and
A. J.
Goddard
, “
Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations
,”
Comput. Methods Appl. Mech. Eng.
190
,
3771
–
3796
(
2001
).
https://doi.org/10.1016/S0045-7825(00)00294-2
Google Scholar
Crossref
Search ADS
 
© 2024 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).
Published open access through an agreement with University of Aberdeen, Saint Mary Kings College