Two-fluid model analyses of
instabilities and non-uniformities in bubbly gas-liquid flows

Henrik Ström^a,*, Klas Jareteg^b,
Srdjan Sasic^a and Christophe Demazière^b

^aDivision
of Fluid Dynamics, Department of Applied Mechanics

^bDivision
of Nuclear Engineering, Department of Applied Physics

Chalmers University of Technology,
SE-412 96 Gothenburg, Sweden

^*Corresponding
author: henrik.strom@chalmers.se

Bubbly gas-liquid flows are of huge importance in
the nuclear industry. The intricate details of the steam-water two-phase flow
in a boiling water reactor significantly affects the coolant properties, which
in turn influences the local heat transfer, the evaporation rates and the transport
of neutrons. Due to the large size of a nuclear reactor, macroscopic models ?
such as the two-fluid model ? are traditionally used to predict the features of
the two-phase flow. Since high-frequency and small-scale phenomena are filtered
out in the derivation of these models, the effect of such phenomena on the
macroscopic fields must be artificially reintroduced, typically via
experimentally-derived closure relationships. It could be speculated that more
general two-fluid models could be developed if the macroscopic formulations
were closed with information from highly resolved simulations with the same
model (much in the vein of the recent development in the area of gas-solid
flows [1-5]). However, if one aspires to abandon the procedure to tune the
macroscopic model directly to the available experimental data, it becomes very
important to first scrutinize the behavior of the basic two-fluid model
formulation at the finest scales. This observation is reinforced by the
knowledge that both the formulation of the two-fluid model and the numerical
algorithms employed to solve it will influence the stability of the solution [6-9].

In the present work, we investigate the performance
of the two-fluid model on smaller scales than what is typically used in the
nuclear industry today. We formulate a simplistic two-fluid model for bubbly
steam-water flow and perform numerical simulations in periodic 2D domains,
investigating the spontaneous emergence of non-uniform bubble volume fraction
fields in the form of meso-scales. We examine how these instabilities are
affected by the formulation of the case specifications, the model formulation
and the methods used to obtain the numerical solution.

The two-fluid model used in the present work is simple
enough to enable direct observations of the effects of each feature added to
the model. It is assumed that the flow can be represented by spherical, rigid
bubbles occupying a certain volume fraction in a continuous liquid. The flow is
isothermal, without any mass transfer, coalescence or breakup. The mass and
momentum balance equations of the two-fluid model employed in the present work
can therefore be formulated as:

Here, i and j are the indices of the
two phases, α is the volume
fraction, r the density, u
the velocity, p the pressure,  the
molecular stress tensor,  the
turbulent stress tensor, g the gravitational acceleration and
K is the momentum exchange coefficient. For the formulation of the
latter, it is assumed that interfacial momentum transfer is entirely dominated
by the drag force, and that this drag force can be characterized by a
correlation developed for rigid spheres. The turbulent stress tensor is either
neglected (assuming that the velocity fluctuations are mainly caused by
large-scale vortical structures [10]) or obtained from the standard k-e
model (assuming that the velocity fluctuations are mainly caused by shear-induced
turbulence [11]).

As a means to quantify the magnitude of meso-scale
structures, we define a global, time-resolved uniformity index such
that:

where  is
the maximum volume fraction of the dispersed phase in the solution domain at
time t,  is
the corresponding minimum volume fraction and  is
the domain-average. The occurrences of non-zero values of  with
time are indicative of the appearance of a non-uniform volume fraction field,
and hence  can
be used as an indicator of the onset of mesoscopic instabilities in an
initially uniform flow.

The temporal development of  is
investigated for a number of computational cases that constitute variations of
a base case (cf. Table 1) inspired by Benyahia and Sundaresan [12]. The domain
is fully periodic, and at time zero, all fields (velocities, pressure and
volume fractions) are uniform and the flow is stationary. The gravitational
acceleration acts in the negative vertical direction and the weight of both
fluids is balanced by a prescribed pressure drop in the same direction.

Table 1. Base case
specification.

A number of investigations of the behavior of  have
been carried out and we here choose to discuss three of them. When varying the average
bubble loading (0.01, 0.05 and 0.1, respectively), a clear trend can be
discerned, as shown in Figure 1. In all cases, the uniformity index starts off
from a very small value, indicative of the uncertainty in the numerical
precision used to store the solution variables. However, as the numerical
simulation proceeds in time, these infinitesimal fluctuations are not dampened
but instead have a tendency to grow. The growth rate of these disturbances
increases until the solution reaches a maximally unstable state. The higher the
volume fraction of the dispersed phase, the sooner this state is attained.
After having reached the maximally unstable state, the solution remains
non-uniform throughout the rest of the simulation, although the magnitude of  decreases
slowly with time.

Figure
1.
The evolution of  for the base case
with three different average bubble loadings.

Physically, the maximally unstable state corresponds
to a striped pattern of the bubble volume fraction field, with alternating
bubble-rich and bubble-lean regions. With the passage of time, these regions are
smoothed out by diffusive mechanisms. Predictions obtained with different
numerical algorithms yield qualitatively similar results. In this context,
?numerical algorithm? represents all possible variations of the procedure
employed to obtain the numerical solution to the posed mathematical problem,
including everything from spatial and temporal resolution to discretization of
the governing equations and the specific algorithms implemented in the
underlying solver. As an example, Figure 2 shows a comparison of two
predictions of the dispersed phase volume fraction field at the maximally
unstable state obtained using different numerical algorithms. The exact magnitude
and appearance of the fluctuations varies but the overall picture is the same.
Similarly, the smoothing out of the striped pattern a few seconds after the
passage through the maximally unstable state produces the same type of behavior
of the time-resolved non-uniformity index, but the details of the meso-scale
structures are not identical. We take this to mean that the exact time-history
obtained is unlikely to reflect a true physical behavior of the system, but
should be interpreted rather as a display of the unstable character of a bubbly
two-phase flow and the complex interplay between the design of the mathematical
model and the choice of numerical algorithms employed in solving it.

Figure
2.
Examples of the dispersed phase volume fraction field at (a) and a few seconds
after (b) the attainment of the maximally unstable state in numerical
simulations with two different numerical algorithms. The average bubble loading
is 0.05.

In many bubbly flows of industrial relevance,
shear-induced turbulence also plays an important role. It is then customary to model
the shear-induced turbulent stress tensor using a Reynolds-Averaged
Navier-Stokes (RANS) turbulence model. The averaging performed in the
derivation of the RANS equations implies that the velocity fluctuations are
filtered out, and it is therefore not surprising that our work has shown that the
effect of including a RANS-based turbulence model is to dampen out non-uniformities.
There is, however, a side effect in that a suppression of the velocity
fluctuations in the two-fluid model by a turbulence model also effectively hinders
the appearance of fluctuations in the volume fraction field by neglecting the correlations
between the two types of fluctuations. The results obtained in the present work
therefore highlight a need to study the details of the interaction between
turbulence and the dispersed phase meso-scale structures using more
comprehensive mathematical frameworks.

ACKNOWLEDGEMENTS

The
research was conducted with funding from the Swedish Research Council
(Vetenskapsrådet) as a part of the Development of Revolutionary and Accurate
Methods for Safety Analyses of Future and Existing Reactors (DREAM4SAFER)
framework grant (contract number C0467701). The Swedish Center for Nuclear
Technology (SKC) is acknowledged for financially supporting K. Jareteg. The
computations were partly performed on resources at Chalmers Centre for
Computational Science and Engineering (C3SE) provided by the Swedish National
Infrastructure for Computing (SNIC).

REFERENCES

[1]
Agrawal et al., J. Fluid Mech. 445 (2001) 151.

[2]
Zhang & VanderHeyden, Powder Tech. 116 (2001) 133.

[3]
Igci et al., AIChE J. 54 (2008) 1431.

[4]
Wang et al., Chem. Eng. Sci. 64 (2009) 622.

[5]
Yang et al., Ind. Eng. Chem. Res.
43 (2004) 5548.

[6] Gudmundsson,
PhD thesis, KTH Royal Institute of Technology, Stockholm, Sweden, 2005.

[7]
Coquel et al., J. Comp. Phys.
136 (1997) 272.

[8] Kreiss & Yström, Math.
Computer Modelling 35 (2002) 1271.

[9]
Sankaranarayanan & Sundaresan, Chem. Eng. Sci. 57 (2002) 3521.

[10]
Ojima et al., Int. J. Multiphase Flow 67 (2014) 111.

[11]
Tomiyama & Shimada, J. Pressure Vessel Tech. 123 (2001) 510.

[12]
Benyahia & Sundaresan, Powder Tech. 220 (2012) 2.