# Two-Fluid Model Analyses of Instabilities and Non-Uniformities in Bubbly Gas-Liquid Flows

text-align:center">
font-family:"Times New Roman","serif"'>Two-fluid model analyses of

instabilities and non-uniformities in bubbly gas-liquid flows

text-align:center"> font-family:"Times New Roman","serif"'>

text-align:center">
font-family:"Times New Roman","serif"'>Henrik StrÃ¶m^{a,*}, Klas Jareteg^{b},

Srdjan Sasic^{a} and Christophe DemaziÃ¨re^{b}

text-align:center">^{
115%;font-family:"Times New Roman","serif"'> }

text-align:center">^{
115%;font-family:"Times New Roman","serif"'>a}*Division
of Fluid Dynamics, Department of Applied Mechanics*

text-align:center">^{
115%;font-family:"Times New Roman","serif"'>b}*Division
of Nuclear Engineering, Department of Applied Physics*

text-align:center">*
115%;font-family:"Times New Roman","serif"'>Chalmers University of Technology,
SE-412 96 Gothenburg, Sweden*

text-align:center">^{
115%;font-family:"Times New Roman","serif"'>*}*Corresponding
author: *

*henrik.strom@chalmers.se*

justify;border:none;padding:0cm"> line-height:115%;font-family:"Times New Roman","serif"'>

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>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].

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>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.

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>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:

justify"> "Times New Roman","serif"'>

justify">

justify"> "Times New Roman","serif"'>

justify">

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>Here, *i *and *j* are the indices of the

two phases, *
115%;font-family:Symbol">α*
line-height:115%;font-family:"Times New Roman","serif"'> is the volume

fraction, *
font-family:Symbol">r*
line-height:115%;font-family:"Times New Roman","serif"'> the density, *u*

the velocity, *p* the pressure,
position:relative;top:5.5pt'>**Â **Â the

molecular stress tensor,
position:relative;top:5.5pt'>**Â **the**
**turbulent stress tensor,

**the gravitational acceleration and**

*g**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]).

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>As a means to quantify the magnitude of meso-scale

structures, we define a global, time-resolved uniformity index
position:relative;top:5.5pt'>such

that:

justify"> "Times New Roman","serif"'>

text-align:center">

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>where
position:relative;top:6.0pt'>Â is

the maximum volume fraction of the dispersed phase in the solution domain at

time *t*,
position:relative;top:6.0pt'>Â is

the corresponding minimum volume fraction and
position:relative;top:7.0pt'>Â is

the domain-average. The occurrences of non-zero values of
position:relative;top:5.5pt'>Â with

time are indicative of the appearance of a non-uniform volume fraction field,

and hence
position:relative;top:5.5pt'>Â can

be used as an indicator of the onset of mesoscopic instabilities in an

initially uniform flow.

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>The temporal development of
position:relative;top:5.5pt'>Â 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.

justify"> "Times New Roman","serif"'>

justify">

justify">**
font-family:"Times New Roman","serif"'>T****able 1.** Base case

specification.

font-family:"Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>A number of investigations of the behavior of
position:relative;top:5.5pt'>Â 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
position:relative;top:5.5pt'>Â decreases

slowly with time.

text-align:center"> "Times New Roman","serif"'>

justify">**Figure
1.**

The evolution of position:relative;top:5.5pt'>Â for the base case

with three different average bubble loadings.

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>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.

justify"> "Times New Roman","serif"'>

text-align:center"> "Times New Roman","serif"'>

justify">**
font-family:"Times New Roman","serif"'> **

justify">**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.

justify"> "Times New Roman","serif"'>

justify">
"Times New Roman","serif"'>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.

justify"> "Times New Roman","serif"'>

justify"> "Times New Roman","serif"'>

**
115%;font-family:"Times New Roman","serif"'> **

justify">**
font-family:"Times New Roman","serif"'>ACKNOWLEDGEMENTS**

justify">**
font-family:"Times New Roman","serif"'> **

justify">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).

justify"> "Times New Roman","serif"'>

justify">**
font-family:"Times New Roman","serif"'> **

justify">**
font-family:"Times New Roman","serif"'>REFERENCES**

justify">**
font-family:"Times New Roman","serif"'> **

[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. **
line-height:115%;font-family:"Times New Roman","serif"'>Chem. Res.*

43 (2004) 5548.

[6] Gudmundsson,

PhD thesis, KTH Royal Institute of Technology, Stockholm, Sweden, 2005.

[7]

Coquel et al., *J. Comp. **
line-height:115%;font-family:"Times New Roman","serif"'>Phys.*

136 (1997) 272.

[8]
line-height:115%;font-family:"Times New Roman","serif"'>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.