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

Authors: F. Moretti; D. Melideo; A. Del Nevo; F. D'Auria; T. Höhne; E. Lisenkov
Journal: Science and Technology of Nuclear Installations
 (2009)

---

## Abstract

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.


---

## Body

## 1. 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
single-phase in-vessel flows. Analogous analyses were performed, for instance,
for some experiments carried out on the above-mentioned 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.


## 2. Description of the Experiment

The experimental facility basically consists of an RPV
model, connected with four circulating loops. The RPV model (Figure1(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.(a) Vertical cross-section of the RPV model; (b) 3D isometric sketch of the facility.
(a)(b)Each loop is equipped with an independent
computer-controlled 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 (Figure1(b)) do not exactly reproduce the real piping layout;
however the related volumes are such that the 1:53 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 (Figure2). Furthermore, a continuous tracer
injection can also be performed into the volume compensation tank located
upstream of a circulation pump.Figure 2
Location of the tracer slug.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 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 m3
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 Q0=220 m3/h);
10 seconds are enough to reach ~98% of the target
flowrate: (1)Q=Q0⋅(1−e−0,25⋅τ). 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.


## 3. Description of the Computational Model

### 3.1. Computational Grid

The computational domain selected for the in-vessel
mixing simulations (shaded region in Figure3) 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.Sketch of the
computational domain chosen for CFD simulations.
(a)(b)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).(viii)
Core support plate.(ix)
Baffle inner wall.
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 high-quality 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 distancey+ 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.Table 1
Size of reference grid (M=106).
No. of nodesNo. of tetrahedraNo. of wedgesNo. of pyramidsNo. of hexahedraTotal no. of elementsReference grid4.3 M3 M630 00032 0002.75 M6.46 MIt 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 Figure4.Reference grid.
(a)(b)(c)(d)


### 3.2. Simulations Setup

The simulations have been performed with the
commercial, multipurpose CFD code ANSYS CFX-10.0 (see [15]), using 8
processors of a Linux-cluster 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/m3,(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 time-steps 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 time-step 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).(ii)
5%
turbulence intensity at loop #4 inlet nozzle.(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.Loop flowrates (posttest values coincide with experimental values).
(a)(b)(c)(d)


## 3.1. Computational Grid

The computational domain selected for the in-vessel
mixing simulations (shaded region in Figure3) 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.Sketch of the
computational domain chosen for CFD simulations.
(a)(b)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).(viii)
Core support plate.(ix)
Baffle inner wall.
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 high-quality 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 distancey+ 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.Table 1
Size of reference grid (M=106).
No. of nodesNo. of tetrahedraNo. of wedgesNo. of pyramidsNo. of hexahedraTotal no. of elementsReference grid4.3 M3 M630 00032 0002.75 M6.46 MIt 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 Figure4.Reference grid.
(a)(b)(c)(d)


## 3.2. Simulations Setup

The simulations have been performed with the
commercial, multipurpose CFD code ANSYS CFX-10.0 (see [15]), using 8
processors of a Linux-cluster 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/m3,(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 time-steps 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 time-step 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).(ii)
5%
turbulence intensity at loop #4 inlet nozzle.(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.Loop flowrates (posttest values coincide with experimental values).
(a)(b)(c)(d)


## 4. 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.Figure6(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 horizontal component and tends to reach the
opposite side before moving downwards (towards the lower plenum). Besides, a
portion of the flow leaves the RPV model through the idle loops (the related
valves being kept open): such inverse flows are shown by some streamlines in
the picture, and are expected to affect the amount of tracer that will reach
the core inlet (since part of the tracer will exit through the idle loops).
Furthermore, a stagnation region appears below the starting loop. This is also
shown by the azimuthal profile of the velocity in the DC (at different instants)
plotted in Figure 6(b).Numerical results: (a) velocity field
(streamlines from loop 4); (b) azimuthal velocity profile in DC.
(a)(b)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.Figure7 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 measurements, 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.Figure 7
Comparison of mixing scalar distribution at
core inlet during slug passage.The pretest results show the same results, from a
qualitative point of view. In particular, the appearance of 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 Figure5, 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 Figure8, 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 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.Figure 8
Maximum mixing scalar at core inlet.Another key parameter is the core-averaged
perturbation, and the related results comparison is shown in Figure9 (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.Figure 9
Core-averaged mixing scalar.A time integration of the core-averaged perturbation
provides a measure of the “accumulated perturbation”; this is shown in Figure10. 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.Figure 10
Accumulated perturbation at core inlet.The accumulated perturbation at 25 seconds for both
the experiment and the posttest results is shown in Figure11 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.Maps of channel-by-channel accumulated
perturbation at 25 seconds.
(a)(b)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 Figure11.Table2 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).Table 2
Comparison of results (perturbation
appearance; max. perturbation; core-average).
Exp−σExpExp+σPretestPosttestTime for MS  >  0.1 (s)—8.5—7.58.7Time of Max. MS (s)—9.9—910.3Max. MS [−]0.6700.7210.7730.7550.758Time of max. Ave. MS (s)—10.9—1011.0Max. Ave. MS (90 ch.) (−)0.3340.3620.3890.4600.489Similar 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 core-averaged 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]): (2)DEV1i,t=ci,t−ei,t, where ci,t and ei,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):(3)DEV2_SIGNi=1tN−t0⋅∑k=1NDEV1i,t⋅(tk−tk−1),DEV2_ABSi=1tN−t0⋅∑k=1N|DEV1i,t|⋅(tk−tk−1),DEV2_RMSi=1tN−t0⋅∑k=1N(DEV1i,t)2⋅(tk−tk−1)2, where N is the number of time-steps within the
selected time period, and tk is the time value at kth time-step.Maps of the DEV2 deviations for both calculations are
plotted in Figure12, 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).Figure 12
Maps of DEV2 deviations.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 Table3 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.Table 3
Core- and time-averaged deviations (DEV3).
PretestPosttestDEV3_SIGN0.00110.0020DEV3_ABS0.06200.0311DEV3_RMS0.00920.0049The local instantaneous deviations can also be
directly averaged over the instrumented locations, so as to obtain
time-dependent deviations (DEV4), according to the following equations:(4)DEV4_SIGNt=1M⋅∑k=1MDEV1i,t,DEV4_ABSt=1M⋅∑k=1M|DEV1i,t|,DEV4_RMSt=1M⋅∑k=1M(DEV1i,t)2. 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.Core-averaged deviations (DEV4): sign, abs.
value, root mean square.
(a)(b)(c)The second and the third plots show the same qualitative behaviour; in both cases the noticeable improvement of posttest results is evident.


## 5. 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.



---

*Source: 436218-2008-12-15.xml*