CFD Analysis of a Slug Mixing Experiment Conducted on a VVER-1000 Model

A commercial CFD code was applied, for validation purposes, to the simulation of a slug mixing experiment carried out at OKB “Gidropress” scaled facility in the framework of EC TACIS project R2.02/02: “Development of safety analysis capabilities for VVER-1000 transients involving spatial variations of coolant properties (temperature or boron concentration) at core inlet.” Such experimental model reproduces a VVER-1000 nuclear reactor and is aimed at investigating the in-vessel mixing phenomena. The addressed experiment involves the start-up of one of the four reactor coolant pumps (the other three remaining idle), and the presence of a tracer slug on the starting loop, which is thus transported to the reactor pressure vessel where it mixes with the clear water. Such conditions may occur in a boron dilution scenario, hence the relevance of the addressed phenomena for nuclear reactor safety. Both a pretest and a posttest CFD simulations of the mentioned experiment were performed, which differ in the definition of the boundary conditions (based either on nominal quantities or on measured quantities, resp.). The numerical results are qualitatively and quantitatively analyzed and compared against the measured data in terms of space and time tracer distribution at the core inlet. The improvement of the results due to the optimization of the boundary conditions is evidenced, and a quantification of the simulation accuracy is proposed.


Introduction
In a pressurized water reactor (PWR) several transient scenarios can be hypothesized leading to a perturbation of the coolant time and space distribution at the core inlet (such as temperature and boron concentration), which in turn can induce positive reactivity insertion and power excursion.Transients leading to Boron dilution as well as main steam line break (MSLB) transients are examples of such scenarios.
The perturbation is influenced by the turbulent mixing phenomena occurring inside the reactor pressure vessel (RPV), that is, a perfect mixing between the perturbed coolant (e.g., a deborated slug coming from a loop) and the nonperturbed coolant is expected to lead to the smallest core response, while the absence of mixing is likely to induce a stronger and localized reactivity insertion.Obviously a quantitative assessment of the relationship between the mixing effects and their consequences in terms of reactivity is needed for demonstrating the reactor safety.
The mixing phenomena are inherently three-dimensional, therefore they can be properly analyzed and predicted by means of numerical tools having 3D capabilities, in particular the computational fluid dynamic (CFD) codes (and, to a certain extent, by system codes embedding 3D modules and mixing models).
Several international projects and experimental campaigns have been conducted in the past to investigate the in-vessel mixing phenomena and the code capabilities to predict them.Examples are the experiments carried out at Forschungszentrum Dresden-Rossendorf (FZD) ROCOM facility [1], University of Maryland [2], Vattenfall [3], while as far as the code assessment is concerned, the OECD/NEA International Standard Problem (ISP) no.43 [4], the EC FLOMIX-R project [5], and the EC ECORA project [6] can be mentioned.
Recently, these issues have been addressed in the framework of the EC-funded TACIS project R2.02/02 "Development of safety analysis capabilities for VVER-1000 transients involving spatial variations of coolant properties (temperature or boron concentration) at core inlet" [7].An extensive experimental campaign was conducted at the OKB "Gidropress" mixing facility (a 1 : 5 scaled model of a VVER-1000 reactor) to study fluid mixing scenarios featured by different flow conditions, such as symmetric and asymmetric steady pump operation at nominal flowrates in presence of tracer injection (5 experiments), and pump start-up scenarios in presence of tracer slugs (5 more experiments).All the measured data collected have been utilized for the validation of mixing models implemented in a set of Russian thermo-hydraulics system codes, with CFD being used as a valuable support to the phenomena understanding, and results interpretation, and being object of validation itself.
One of the experimental tests performed consisted in a main coolant pump (MCP) start-up, with the other three pumps switched-off and a tracer slug (simulating deborated water) accumulated in the cold leg of the starting loop.Such experiment was then simulated both with system codes and CFD codes.In particular, pretest and posttest simulations were run using the commercial CFD code ANSYS-CFX and the results obtained were compared against the measured data.Such CFD code validation activity is described in the present paper.
As explained above, the present work is a part of a wider and more comprehensive activity, which included the CFD grid generation, the pretest and posttest simulations of all the experiment performed, the execution of sensitivity analyses on the main modelling parameters, in compliance with the requirements of the Best Practice Guidelines (BPG, [8,9]).Description of the entire work is beyond the scope of the paper and is not reported, but some additional information can be found for instance in [10,11].
This work is connected to the CFD code validation activity in progress at the University of Pisa, related to singlephase in-vessel flows.Analogous analyses were performed, for instance, for some experiments carried out on the abovementioned ROCOM facility [12].
It is also worth mentioning that CFD validation activities [13] had been carried out in the recent past on a previous version of the same Gidropress mixing facility in the framework of the above-mentioned FLOMIX-R project.

Description of the Experiment
The experimental facility basically consists of an RPV model, connected with four circulating loops.The RPV model (Figure 1(a)) is made of steel and reproduces, at a 1 : 5 scale, practically all the geometrical features of the RPV of a VVER-1000 reactor (namely, Novovoronezh NPP reactor, Unit no. 5) which are affecting the in-vessel mixing phenomena up to the core inlet, particularly the internal components such as the barrel, the lower ellipsoidal perforated shell (with more than 1300 drillings of two different diameters), the core support columns (one for each of the 151 fuel assemblies) and the core lower plate (which separates the core region from the lower plenum region).The core region is actually not modelled; rather a structure is present made of perforated plates and guide tubes supporting 90 conductivity probes which are located just above the lower core plate.Each loop is equipped with an independent computercontrolled circulation pump, which permits to simulate a wide range of flow conditions.An expansion tank is connected to one of the loops in lieu of the pressurizer; atmospheric conditions reign above the water level in the tank, while higher pressures (although still in the range 1 to 2 atm) occur in the RPV model due to the hydrostatic effect and to the pumps head.The experiments are conducted at ambient temperature.The circulating loops (Figure 1(b)) do not exactly reproduce the real piping layout; however the related volumes are such that the 1 : 5 3 volume scale is kept.Some auxiliary systems are present for the tracer injection, consisting in injection pumps, fast acting valves, a tracer tank, and pipelines connecting all such components to the main loops.For instance, such systems can be operated for accumulating a tracer slug in the ascending section of one loop while the pumps are at rest, with such section being "isolated" by two fast acting valves (Figure 2).Furthermore, a continuous tracer injection can also be performed into the volume compensation tank located upstream of a circulation pump.
The tracer utilized is sodium chloride, which alters the water electrical conductivity.Through a calibration procedure, the conductivity can be easily correlated to the salt concentration.The facility is equipped with a number of conductivity probes, providing high-frequency measurements of the local tracer concentration.As mentioned above, 90 of such probes are located above the lower core plate, each being aligned with the centreline of one coolant channel.This means that experimental information is available for 60% of the coolant channels (90 out of 151), which is obviously not an "ideal" configuration (as it would be if all the Science and Technology of Nuclear Installations channels were instrumented); however such measurement equipment still permits to gather valuable information of the perturbation at the core inlet.A conductivity probe is located also at each inlet and outlet nozzle; some probes are present in the tracer tank.
The loop flowrates are measured by electromagnetic flow meters located close to each inlet nozzle.
As can be understood from the description above, the facility can be operated such as to simulate a wide spectrum of operation conditions and accidental scenarios involving the perturbation of the coolant properties distribution at the core inlet.The experiment addressed in the present work was intended to reproduce the start-up of one reactor coolant pump (the other pumps remaining at rest) assuming that a "deborated slug" had previously been accumulated in the starting loop.The slug is thus transported inside the RPV, where it partially mixes with the normally borated water before reaching the core inlet and then introducing a positive reactivity in the reactor core.
Namely, the deborated slug is here represented by a salted water slug (0.072 m 3 volume, which roughly corresponds to the scaled volume of the loop seal, where a deborated slug would most probably accumulate).
The starting pump is run, via the numerical control, such as to achieve an exponential growth for the flowrate, according to (1) (the target flowrate being Q 0 = 220 m 3 /h); 10 seconds are enough to reach ∼98% of the target flowrate: The isolation valves of the idle loops are left open; therefore inverse flows develop which are expected to strongly affect the flow field in the RPV model as well as the tracer distribution.The inverse flowrates are not known before the execution of the experiments, and thus constitute the main unknown parameters in the pretest phase of the numerical analysis.

Description of the Computational Model
3.1.Computational Grid.The computational domain selected for the in-vessel mixing simulations (shaded region in Figure 3) includes the following coolant regions: cold legs, inlet nozzles, downcomer (DC), lower plenum (LP).The reactor core region and the upper plenum are not modelled because they are not expected to influence the coolant flow upstream of the core inlet.However, a dummy outlet volume is defined corresponding to a fraction of the core region, to permit the easy application of pressure-controlled outlet boundary conditions.The identified computational domain is defined and bounded by the following solid parts.
(i) Inner wall of the cold legs and inlet nozzles (including round surface at connection with vessel wall).
(ii) Inner wall of the vessel (including cylindrical regions, diameter variations, elliptical bottom).
(iii) Consoles, located in the lower part of the DC.
(iv) Outer wall of the barrel (including elliptical bottom).
(v) Inner wall of the barrel (including elliptical bottom), only up to the core inlet.
(vi) Holes through the barrel bottom (also referred to as "perforated shell" in the following).
(vii) Support columns, located in the region between the inner wall of the barrel bottom and the lower side of the core support plate; each column includes a "solid column" part (14 mm diameter) on the bottom and a "perforated column" part on the top (a tube, 38 mm outer diameter, connected to the solid columns through a conic region, and having perforations on its wall allowing the fluid to pass from the LP to the core support plate holes and then to the core region).
The presence of such a large number of small geometric details (consoles, perforations through the barrel bottom, support columns, etc.) makes the achievement of a highquality and accurate computational grid quite a tough task.
The mesh has been developed with the package ANSYS ICEM-CFD 10.0 (see [14]), following a modular approach, that is, the domain has been subdivided into several subdomains which have been meshed separately.Then the submeshes obtained have been connected together by means of "interfaces" (one-to-one interfaces were used for conformal subgrids and general grid interfaces for nonconformal subgrids; see [15]).This approach allowed adopting different mesh types in different subdomains, namely the DC was meshed with hexahedral elements, while tetrahedrons were used to model LP region (where the complexity of the geometry would make impracticable the hexahedral meshing), including the ellipsoidal perforated shell with its ∼1300 holes.The result is a so-called hybrid grid.
For a proper treatment of the near-wall turbulence, based on logarithmic wall functions, the grid spacing was refined close to the walls in the hexahedral submeshes, while prism layers were inflated in the tetrahedral submeshes.The adopted spacing yielded an average value of about 110 for the nondimensional distance y+ of the first nodes from the walls.
Several grids were generated and assembled based on different meshing approaches (e.g., tetrahedral versus hexahedral elements) and sizes adopted in some submeshes.Grid sensitivity analyses are described in [11].They helped selecting a reference grid (the same used for the present calculations; see Table 1) as the one providing the better convergence, and show that an improvement of the results could be expected from finer meshes.
It is worth remarking that the reference grid is to be considered as a "production grid," in the sense that its size results from a compromise between the need of achieving a high numerical accuracy and mesh-converged results (as recommended by the BPG) on one side, and the computational resources limitations on the other side.It was not possible to demonstrate that the grid is able to provide grid-independent results, as usually happens when addressing CFD problems having the same degree of complexity.However, it is believed to be a state-of-the-art grid, suitable for CFD simulation of turbulent flows, at least as far as the Reynolds-Averaged Navier-Stokes turbulence modelling is adopted.Some pictures of the reference grid are shown in Figure 4.

Simulations Setup.
The simulations have been performed with the commercial, multipurpose CFD code ANSYS CFX-10.0 (see [15]), using 8 processors of a Linuxcluster available at the University of Pisa.The main features of the simulations setup are as follows: (i) working fluid: water (incompressible) at 1 atm, 25 • C, (ii) density: 997 kg/m 3 , (iii) dynamic viscosity: 8.899 × 10 −4 kg m −1 s −1 , (iv) turbulence accounted for with SST model.
The following field equations have been solved: (i) mass balance (continuity), (ii) momentum balance (Navier-Stokes), (iii) transport of turbulent kinetic energy (k), (iv) transport of turbulent eddy frequency, (v) transport of an additional, user-defined, scalar variable simulating the tracer.The tracer concentration is handled in terms of normalized concentration (also referred to as the mixing scalar or MS).Normalization is such that the mixing scalar ranges between the values 0 and 1, which correspond respectively to absence of tracer (i.e., full boron concentration in a hypothetical real plant transient) and initial concentration in the tracer slug (i.e., lowest boron concentration).
Since the addressed flow is dominated by turbulent diffusion, the molecular diffusion of the tracer provides a negligible contribution to the effective diffusion and was then neglected.Sensitivity analyses on the tracer diffusivity performed within the FLOMIX-R project (see [5]) support this assumption.
The transient solver available in CFX was used for both calculations, and the second-order backward Euler time advancement scheme was adopted.A constant time-step equal to 0.05 second was used, which made most timesteps converge with only two internal iterations, and allowed obtaining the results in a reasonable time (about 10 days computation on 8 CPUs of a AMD Opteron Linux cluster, for simulating 25 seconds).Sensitivity analyses on the timestep size are envisaged for the future.
The upwind scheme for the discretization of the advection terms was selected; adopting higher-order schemes is generally recommended (see, e.g., the Best Practice Guidelines [8]), because they are less prone to numerical diffusion than first-order schemes (such as upwind), however previous sensitivity calculations performed using the same grid had shown some nonsatisfactory performance (local nonphysical oscillations, bad convergence) when a higher-order scheme was used, therefore it was decided to stay with the upwind scheme.
The initial conditions (for both the pretest and the posttest calculations) consisted in zero-velocity flow over the whole domain, and zero-concentration everywhere except for the volume corresponding to the tracer slug, which was marked with mixing scalar equal to 1.The following boundary conditions were set for the pretest calculation.(i) Time-dependent flowrate at loop #4 inlet nozzle, according the theoretical law (see (1) above).
(iii) Pressure-controlled "Opening" at inlet nozzles #1, 2 and 3 (to permit inverse flows), with additional concentrated pressure losses to account for the overall flow resistance of the idle loops (the pressure loss coefficients have been roughly estimated based on sensitivity calculations and experimental information on the inverse flowrates, which were known not to exceed 10% of the nominal flowrate).(iv) Pressure-controlled "Outlet" at the top boundary of the dummy outlet volume replacing the core region.
(v) No-slip condition at all walls (i.e., all boundaries not mentioned above).
(vi) Near-wall treatment of turbulence based on logarithmic law.
The posttest calculation setup is identical to the pretest, except for the boundary conditions at the cold legs.In this case, in fact, all the flowrates (including the inverse ones) were imposed based on the measured values.The experimental flowrates are plotted in Figure 5, along with those resulting from the pretest calculation.

Results
All results and experimental data are reported (and compared) in terms of normalized concentration (mixing scalar) at the core inlet, in particular at the 90 instrumented channels locations.Figure 6(a) provides a picture of the flow pattern developing in the RPV model, by means of streamlines entering from the starting loop.The entering flow keeps a dominant Such a qualitative behaviour of the flow is highly dominated by three-dimensional features, so that it would be hardly described by system codes (even if with 3D capabilities).CFD codes represent the "natural" approach to deal with such behaviour, although an accurate modelling of the turbulence may still be a challenging task due to the high anisotropy of the turbulence parameters expected in a strongly bounded flow.
The correct description of the flow field developing in the downcomer is important because it determines the space distribution of the perturbation at the core inlet, particularly the location and the shape of the perturbation.
Figure 7 shows a qualitative code-to-experiment comparison of the mixing scalar at the core inlet at several selected instants during the slug passage.As from the experimental the first perturbation appears at the core inlet at around 9 seconds and is located below loop no. 1, that is, on the opposite side to the starting loop (i.e., no.4).Then the perturbation extends to other peripheral channels in the clockwise direction; furthermore, a secondary perturbation spot appears just below loop no. 2. After a couple of seconds from the first perturbation appearance, almost all channels are affected, and the mixing scalar distribution has become relatively uniform.In a few more seconds, the perturbation disappears from the core inlet.
The pretest results show the same results, from a qualitative point of view.In particular, the appearance a primary perturbation on the opposite side with respect to the starting loop and a secondary perturbation spot below the same loop is correctly described, although with a small discrepancy in timing (1 second ahead) and somewhat larger spatial gradients.Moreover, when most of the perturbation is crossing the core inlet, the spatial distribution is quite less uniform than observed in the experiment.
As can be observed in Figure 5, the pretest calculation overestimated all the inverse flowrates in idle loops.In addition, also the direct flowrate in the starting loop is larger that the measured value (as the experiment is not exactly behaving according to the theoretical law), and this explains why the perturbation reaches the core inlet in advance with respect to the test.Such time shift disappears in the posttest calculation, where the experimental loop flowrates are imposed as boundary conditions.The perturbations appearance now appears aligned with the experiment.
It is evident how the morphology of the perturbation affecting the core inlet is determined by the flow distribution in the downcomer (described above).
It is also evident that the predicted spatial distribution of the perturbation is quite less uniform than observed in the experiment.In other words, a less effective mixing is predicted, as it has previously been observed in similar works (see [12]), and this behaviour is most probably related to limitations of the RANS turbulence modelling.
A key parameter affecting the core neutron kinetics response is the maximum perturbation (e.g., the lowest boron concentration, in a boron dilution scenario) reached at the core inlet.The related code predictions are plotted in Figure 8, where they are compared with the corresponding experimental trend.As mentioned before, five runs were conducted for this experiment, and the measured values were averaged over such data sets.The mean value of the maximum perturbation is reported in the figure, along with the two curves defining a confidence interval of one standard deviation around the mean value.It is observed that both calculations slightly overpredicted the peak of the mean value curve, although still within the confidence interval.The pretest results show a time shift of 1 second in advance (as already observed from Figure 7), which is related to Science and Technology of Nuclear Installations the nonoptimized boundary conditions.The posttest results show instead an accurate timing for the peak occurrence, as well as for the first appearance of the perturbation (around 8 seconds).Later, the code prediction shows a slight delay in the maximum perturbation decrease: at 15 seconds the predicted value for the maximum perturbation is around 0.2, while the experimental value is a little above 0.1.The posttest results, although they are generally outside the confidence interval, look pretty close to the experimental behaviour.
Another key parameter is the core-averaged perturbation, and the related results comparison is shown in Figure 9 (the averaging is made on the 90 instrumented locations, both for measured and calculated data).The pretest results show the same time shift observed above.The posttest results show a correct timing, and a less smooth behaviour than the experimental trend, which indicates that a less "diffused" slug is passing through the core inlet.
A time integration of the core-averaged perturbation provides a measure of the "accumulated perturbation"; this is shown in Figure 10.Again, the posttest results show a less diffusive trend (indicated by a steeper gradient); in other words, the most amount of perturbation takes-according to the code prediction-a smaller time to cross the core inlet than in the experiment.Quite surprisingly, at the end of the slug passage both calculations predicted the same accumulated perturbation as the experiment, thta is, the same amount of tracer has reached the core inlet despite the nonaccurate boundary conditions in the pretest.
The accumulated perturbation at 25 seconds for both the experiment and the posttest results is shown in Figure 11 for each instrumented channel.Those maps evidence that the code tends to underpredict the overall perturbation in the central region and to overpredict it in the peripheral region between loops I and IV.
The measured maximum local accumulated perturbation is 2.02 (s/−), while the predicted value is 2.07 (s/−).The locations of those two maxima are indicated by red circles in Figure 11.
Table 2 summarizes the results obtained for some key parameters such as the timing of perturbation appearance (defined as MS = 0.1), timing and value of the maximum perturbation, and timing and value of the core-averaged perturbation peak.As observed before, the appearance of the perturbation is predicted 1 second in advance by the pretest calculation, and with a 0.2 second delay by the posttest calculation.Similar time discrepancies (−0.9 second and +0.4 second, resp.)appear for the prediction of maximum.The maximum value is predicted quite satisfactorily in both cases (with a 5% overestimation, which is, however, within the ±σ confidence interval).
Similar time discrepancies (−0.9 second and +0.1 second, resp.) also appear for the prediction of core-averaged peak, while the related peak value is noticeably overpredicted in both cases (27% and 35%, resp.).This seems to indicate a less effective mixing.
A quantitative analysis of the agreement between code predictions and measured data requires taking into account the results channel by channel, in addition to the coreaveraged and maximum perturbations discussed above.However, as easily expected from the quantitative analysis shown before, an excellent agreement would be observed at some locations while at the other locations the perturbation will be either overpredicted or underpredicted by the calculations.This does not allow an easy judgement on the overall quality of the code prediction, unless some general, synthetic accuracy parameter is defined.
A local instantaneous code-to-experiment deviation can be defined as follows (based on the same approach adopted within the FLOMIX-R project [1]): where c i,t and e i,t , respectively, represent the calculated and experimental values at ith location and tth time-step.The deviation DEV1 can be averaged over a time interval of interest (e.g., 0-17 seconds, corresponding to the slug passage through the core inlet).The following three deviations are thus obtained (based, resp., on relative and absolute values of DEV1 deviations, and on a root mean square averaging approach): where N is the number of time-steps within the selected time period, and t k is the time value at kth time-step.Maps of the DEV2 deviations for both calculations are plotted in Figure 12, obviously for the 90 instrumented channels only (the others being represented by white colour).Concerning the deviations with their sign, they approximately range between −0.04 and 0.04, and no evident change is observed from pre-to posttest: this is because the two calculations actually behave similarly, except for the time shift, and thus errors with opposite sign during the transient partly compensate.Some locations are evidenced, in both cases, where the perturbation is systematically overpredicted (red) or underpredicted (blue).
Concerning the absolute deviations, it is not possible to identify specific patterns on the map of pretest results, while on posttest map it is observed that the largest discrepancies occur in the central region and in the peripheral region around 90 • away from loop #4 (on both directions); moreover, a noticeable improvement is noticed from pretest to posttest.The same behaviour is observed for the root mean square deviations.
If the DEV2 deviations are averaged over the instrumented locations, then the results in Table 3 are obtained (deviations DEV3), which represent a measure of the overall accumulated deviations.Again, the higher accuracy of the posttest predictions is evidenced.Only the DEV3SIGN deviation is increased.
The local instantaneous deviations can also be directly averaged over the instrumented locations, so as to obtain The resulting plots are shown in Figure 13.The first plot clearly indicates that the pretest results first over predict the perturbation (until 11 seconds), then the under prediction prevails; this is related to the time shift.The posttest results show an opposite behaviour, and generally the discrepancy is much reduced.The second and the third plots show the same qualitative behaviour; in both cases the noticeable improvement of posttest results is evident.

Conclusions
A pump start-up experiment with the presence of a tracer slug, conducted on a Gidropress mixing facility in the framework of TACIS project R2.02/02, was simulated with the CFD code ANSYS CFX.Both a pretest and a posttest calculation were run, differing by the boundary conditions imposed in terms of loop flowrates.The numerical results were compared against the experimental data available, which consist in tracer concentration measurements at several locations at the core inlet.
The results of both calculations showed quite a good agreement with the experiment from the qualitative point of view: in particular, the morphology of the tracer concentration distribution at the core inlet was correctly described, including the appearance of two different perturbation patterns (one on the opposite side with respect to the starting loop, and a secondary one on the same side).The only noticeable difference between the pretest and the posttest-confirmed also by the quantitative analysis-is a time shift (in advance) of the former, due to an imposed loop flowrate which was little higher than actually obtained in the experiment.This qualitative agreement is quite an important achievement, since the addressed scenario is featured by a complex, highly three-dimensional, flow distribution in the downcomer, and its accurate numerical prediction is not a trivial task, due to the well-known limitations of the turbulence modelling based on the Reynolds-Averaged Navier-Stokes approach and particularly on the eddy viscosity concept (i.e., the difficulties in dealing with turbulence anisotropy, unsteady flows, separation phenomena, secondary motions), and typical of most industrial-scale CFD applications.
From a quantitative point of view, the results in terms of maximum perturbation (and related timing), core-averaged perturbation, and accumulated perturbation are also satisfactory.The perturbation peak is overpredicted by 5%, which is comparable with the experimental uncertainty.The predicted time history of the core-averaged perturbation shows a less smooth trend than the experiment, which seems to indicate a less effective mixing (this would be consistent with results from previous CFD validation studies against symmetric loop operation experiments, which had shown a tendency to underpredict the turbulent mixing by the CFD/2-equation turbulence modelling approach).
A further quantitative analysis of the results was done based on a set of "deviations" defined according to a similar approach to that adopted within the FLOMIX-R project.This kind of analysis of the agreement between code predictions and experiment provides a valuable tool to compare the accuracy of different code results.However, a real judgement on the results accuracy cannot be given because it would require a sort of "acceptance thresholds" (in relation to the nuclear reactor safety), which however have not been proposed yet.This is certainly an important matter for future research.
Possible future developments of the present work involve developing finer grids (as far as allowed by the available computing resources), running further sensitivity analyses (e.g., with respect to time discretization, wall roughness) and switching to large eddy simulation (LES) or LES/RANS hybrid approaches for a more accurate prediction of turbulence.

Figure 1 :
Figure 1: (a) Vertical cross-section of the RPV model; (b) 3D isometric sketch of the facility.

Figure 2 :
Figure 2: Location of the tracer slug.

Figure 3 :
Figure 3: Sketch of the computational domain chosen for CFD simulations.