.tifpng.pngconvert #1 .tif

# Large Scale GPU Based Simulations of Turbulent Bubbly Flow in a Square Duct

Purushotam Kumar<EMAIL_ADDRESS>Surya P. Vanka<EMAIL_ADDRESS>
###### Abstract

Turbulent bubbly flows in non-circular ducts are applicable in several
engineering applications, such as nuclear reactors, air conditioning and
refrigeration equipment, thermal power units, electronics cooling etc. In this
paper, we present the results of a numerical study of air-water turbulent
bubbly flow in a periodic vertical square duct. The study is conducted using a
novel numerical technique which leverages Volume of Fluid method for interface
capturing and Sharp Surface Force method for accurate representation of the
surface tension forces. A three-dimensional geometry construction method is
employed during solution of interface equation which gives absolute
conservation of mass and sharp interface between gas and liquid phases. The
spurious velocity which is ubiquitous in numerical modeling of multiphase
flows is substantially smaller with our technique. The entire algorithm has
been implemented on a data parallel mode on multiple graphics processing units
(GPU) taking advantage of the large number of available cores.

We have studied the dynamics of a swarm of spherical bubbles co-flowing with
the upward turbulent flow and compared results with an unladen turbulent flow.
The frictional Reynolds number of the unladen $Re_{\tau}$ is 360, which is
sufficient to sustain a turbulent flow. We observe the turbulence-driven
secondary flows in the mean flow, with complex instantaneous turbulent
vortical structures. The interaction of these secondary flows with the upwards
rising bubbles is very complex and leads to significant changes in the
instantaneous and time-averaged flow field. In this study, we have considered
864 bubbles rising in the square duct on a grid consisting of 32 million
control volumes. The calculations have been for over 1,500,000 time-steps and
took approximately 6 months of calendar time on a 4-GPU computer.

We present the results of mean void fraction distribution, mean velocities,
longitudinal and transverse turbulence intensities along the wall, corner
bisector, and wall bisector. A peak in the void fraction distribution near the
walls is observed representing the migration of bubbles to a preferred section
of the duct. The effects of turbulence-driven secondary flows and
instantaneous large eddies on preferential concentration of the bubbles are
discussed. The dispersed bubbles are seen to break the long elongated
turbulent structures commonly observed in the unladen turbulent flow.

Keywords: Gas-liquid flows; bubble induced turbulence, bubbly flow; VOF;
secondary flow, square duct, GPU computing

##  1 Introduction

The bubbly flow in a circular pipe or a square duct is one of the simplest
fully-developed two-phase flow to envision, yet it provides important
information related to the interaction between bubbles and the continuous
phase. As a consequence, a significant number of experiments have been
conducted to understand upward and downward bubbly flows in circular pipes. In
these studies, the void fraction distribution, mean phase velocities, mean
bubble size, turbulence quantities, wall shear, heat transfer from wall to
fluid etc. have been explored. Many physical quantities of interest such as
heat transfer from wall to fluid, pressure drop and local void fraction depend
on the distribution of bubbles in the duct. Through experiments and numerical
simulations, it has been observed that the bubbles move in the lateral
direction to either accumulate near the walls, also known as wall-peaking
case, or in the central region, referred to as core-peaking case. Intermediate
void fraction distributions with multiple peaks along the radius are also
reported in the literature. This lateral migration of bubbles is related to
the lift force due to fluid-shear and mutual interaction between bubbles,
which are dependent on the bubble size and drift velocity of the dispersed
phase. In the upward bubbly flow, smaller bubbles $(d\leq 3)$ mm concentrated
in the wall region, whereas larger bubbles $(3<d<7)$ mm migrated to the pipe
center. The change in preferential position of bubble is related to the
reversal of lift force [1] as the bubble size increases. Bubbles larger than 7
mm, initially migrate to the central region, but due to shape deformations and
volumetric expansion break-up into smaller bubbles of varying sizes, and an
intermediate void fraction distribution is observed. The lift direction also
changes due to change in the sign of drift velocity. For example, in the
downward bubbly flow, bubbles with diameter smaller than 3 mm migrate to pipe
center, which is in contrast with the upward bubbly flow scenario. The works
of Afghan [2], Zun et al. [3], Liu [4], Liu and Bankoff [5, 6], Tomiyama et
al. [1], Colin et al. [7] can be consulted for further information on void
fraction distributions.
The turbulence characteristics of the surrounding flow can also plays
significant role migration of the bubbles. However, apart from few
experimental and numerical studies to characterize turbulence in channel and
pipe flows, there is very little known about the flow field and turbulence
characteristics of bubbly flow in the square ducts. In general, both mean and
secondary flows are affected by bubbles in the turbulent bubbly flow. This is
related to the preferential distribution of the bubble and subsequent
variation of density in the cross-section. As mentioned earlier, bubbles have
a tendency to concentrate either near the wall or in the central region,
hence, the axial velocity depends strongly on flow conditions. For example, in
the wall peaking scenario the average density near the wall decreases,
therefore for the same pressure gradient the velocity near the wall increases.
This results in a smoothed ‘M’ shaped velocity profile with a flatter central
portion. Similarly, for the core peaking scenario the average density in the
core region decreases, hence, the mean velocity in the central region
increases.
The turbulence generated by bubbles can also have significant impact on the
structures of surrounding fluid flow, and source of this turbulence is
combination of the relative motion between bubbles and fluid, their
deformation, and vortex shedding behind the bubbles. Depending on the bubble
Reynolds number, these vortices can be attached or detached, and their
interactions with other bubbles create fluctuations that are smaller in scale
than that of wall-induced turbulence. Lance and Bataille [8] studied the
interaction between the grid-generated isotropic turbulence and the wall-
induced turbulence in a uniform bubbly flow. They showed that the turbulence
remains isotropic except near the injection point of bubbles. There are two
regimes of interaction between the wall-induced and bubble induced turbulence,
and it depends strongly on the void fraction $\bar{\alpha}$ and
$u_{0}^{\prime}/U_{R}$, the ratio between single phase longitudinal
fluctuations and mean drift velocity. The excess turbulence kinetic energy due
to bubbles first increased linearly with the void fraction for all grid-
generated turbulence intensity, before increasing with a higher slope for
higher grid-generated turbulence intensity. The transition between two regimes
depends on the critical void fraction $\bar{\alpha}\leq 1\%$; below which the
hydrodynamic interactions between the bubbles were negligible, and the
turbulent kinetic energy of the liquid was a summation of the grid-generated
and bubble-induced turbulent kinetic energies. Above the critical void
fraction, the bubbles transferred a great amount of energy to the liquid
phase, primarily due to the mutual interaction between bubbles. Since, the
average distance between the bubbles decreases with the increase in the void
fraction, the hydrodynamic interaction is also expected to increase.
Similarly, when the grid-generated turbulence intensity is increased large
trajectory fluctuations and an increase in probability of mutual interaction
between bubbles is observed. This is probably the reason why excess kinetic
energy increases with the increase in the void fraction and turbulence
intensity of the dispersed phase.
The void fraction also has a strong influence on the longitudinal velocity
fluctuation correlations and the one-dimensional energy spectrum. The
turbulence kinetic energy in the single phase flow is concentrated in low
frequency eddies, and high frequency is described by the typical
$-\frac{5}{3}$ power-law equation. The introduction of bubbles modifies the
turbulence correlation and energy spectrum which can be divided in two
regimes: (a) bubble-induced, and (b) wall-induced. In the first regime, the
shape of correlation and energy spectrum is independent of the void fraction,
but a minor deviation from single phase behavior is observed. In the second
regime, the classical $-\frac{5}{3}$ power-law description of the energy
spectrum is swiftly replaced by $-\frac{8}{3}$ power-law as the void fraction
is increases for at a constant continuous phase flow rate. Hence, the
turbulent kinetic energy in the bubbly flow is spread over a wider range of
frequencies. Using the asymptotic expansion of the turbulence dissipation and
production due to particles, Yuan and Michaelides [9] showed that smaller
particles, whose interaction time $\tau$ is much larger than particle response
time $\tau_{p}$, attenuated the continuous phase turbulence and the reduction
was proportional to $d_{p}^{3}$. On the other hand, the larger particles
($\tau<<\tau_{p}$) enhanced the turbulence which is related to the vortices
generated by the particles. The shedding of vortices were only observed at
particle Reynolds number ($Re_{p}$) larger than 600. The enhancement of
turbulence due to particles was observed at a lower $Re_{p}$ which affirms
role of attached vortices in turbulence enhancement. For a dilute bubbly flow
($\bar{\alpha}\sim 10^{-4}$), Colin et al. [7] showed that the turbulence
intensities enhance in the upward flow and attenuate in the downward flow. In
the upward flow, smaller scales possess more energy in the core region where
as the larger scales possess more energy in the wall region. In the downward
flow, transfer of energy from large to small scales were negligible, and
bubbles always dissipated energy from the liquid.
Hosokawa and Tomiyama [10] studied bubble-induced turbulence in a laminar pipe
flow (bulk Reynolds number = 900) and showed that turbulence kinetic energy,
due to fluctuations introduced by bubbles, increased with the increase in the
void fraction. For a core-peaking case, the turbulence intensities decreased
with distance away from center. For a transition case, the peaks were observed
at $r/R=0.6$, (roughly speaking the trend of fluctuations followed the void
fraction distribution). The radial and circumferential fluctuations were
similar in magnitude and almost half of that in the axial direction. Alméras
et al. [11] studied the bubble-induced turbulence in a large square duct with
void fraction ranging from 1% to 13%. The turbulence intensity in the
streamwise direction was larger than those in the lateral directions. It is
related to the wake of bubbles as they preferentially promote the capture and
transport of fluids in the streamwise direction than in the lateral
directions. They injected a low-diffusive passive dye to elaborate the capture
and transport process. They observed that in both streamwise and lateral
directions, there are two regimes of diffusion coefficient dependence on the
void fraction. At lower void fractions, where hydrodynamic interactions
between bubbles are small, the diffusion coefficient increased nearly as
$\bar{\alpha}^{0.4}$, whereas at higher void fractions, the diffusion
coefficient were seen to be independent of $\bar{\alpha}$. The transition
between two regimes occurred earlier in the lateral directions
$(1\%\leq\bar{\alpha}a\leq 3\%)$ than in the streamwise direction
$(3\%\leq\bar{\alpha}\leq 6\%)$.
In the present paper, we have studied modification of turbulence in a vertical
square duct in upward turbulent bubbly flow. Despite a reasonably large number
of previous studies, the bubble induced turbulence has not been fully
characterized in low void fraction regimes. In particular, there are very few
studies that have examined turbulence modification in square duct. A volume of
fluid method with a sharp surface force method for the inclusion of surface
tension forces has been developed for a multiple graphics processing units
platform. The mean quantities and turbulence statistics are investigated and
comparison between unladen and bubble laden turbulent flow in presented. In
section 2 the governing equations and solution procedure used for the
computations are presented. The computational details and problem description
are presented in Section 4. In Section 5, we present results of the
simulations performed in the study and discuss the important findings. Section
6 presents a summary of the present results.

##  2 Numerical method

We have developed a numerical procedure that solves the relevant governing
equations on a fixed Eulerian collocated grid. The procedure consists of
algorithms for interface capturing, inclusion of surface tension forces, and
solution of Navier-Stokes equations. The volume of fluid (VOF) method of Hirt
and Nichols [12] is used for interface capturing. The truncated interface
between the gas and liquid in a control volume is represented as a oblique
plane, and the interface normal and curvature are computed using a smoothed
liquid volume fraction and the height function method (Rudman [13] and Cummins
et al. [14]), respectively. We have incorporated the surface tension force in
Navier-Stokes equations using a sharp surface force (SSF) and pressure balance
methods (Francois et al.[15], Wang and Tong [16]). The Navier-Stokes equations
are solved using the fractional step method. The numerical procedure has been
previously used to study the effects of duct confinement on dynamics of bubble
rising in a square duct [17, 18], the effects of magnetic field [19] and
variable viscosity [20, 21] on bubble dynamics.

###  2.1 Governing equations

We assume that both the gas and liquid are isothermal, incompressible and
Newtonian fluids. A single fluid approach with a method to capture the
interface between gas and liquid is used. The combined governing equations for
both fluids are given by following equations::

$\displaystyle\nabla\cdot\mathbf{u}=0$ (1)
$\displaystyle\frac{\partial\left(\rho\mathbf{u}\right)}{\partial
t}+\nabla\cdot\left(\rho\mathbf{uu}\right)=-\nabla p$
$\displaystyle+\nabla\cdot\left(\mu\left[\nabla\mathbf{u}+\nabla\mathbf{u}^{T}\right]\right)$
(2)
$\displaystyle+\rho\textbf{g}+\sigma\kappa\mathbf{n}\delta\left(\mathbf{x}-\mathbf{x}_{f}\right)$

In the above equations, $\mathbf{u}$ is fluid velocity, $p$ is pressure with
hydrostatic part being subtracted, $\rho$ is density, $\mu$ is dynamic
viscosity, $\sigma$ is surface tension coefficient, $\kappa$ is interface
curvature, $\mathbf{n}$ is interface normal, $\delta$ is the delta function,
$\mathbf{x}$ is the spatial location where the equation is solved,
$\mathbf{x}_{f}$ is the position of the interface and $\mathbf{g}$ is
acceleration due to gravity.

Since free surface is captured by the volume of fluid method, a time evolution
of VOF function (liquid volume fraction, $\alpha$) given by

$\displaystyle\frac{\partial\alpha}{\partial t}+\textbf{u}\cdot\nabla\alpha=0$
(3)

is solved. The geometry construction method (Rider and Kothe [22]) coupled
with a second-order operator split method (Noh and Woodward [23], Li [24],
Ashgriz and Poo [25] and Francois et al. [15]) is used to solve the evolution
equation. The mixture density and viscosity are calculated as linear weighting
of the individual phase values, as

$\displaystyle\rho$ $\displaystyle=\alpha\rho_{l}+(1-\alpha)\rho_{g}$ (4)
$\displaystyle\mu$ $\displaystyle=\alpha\mu_{l}+(1-\alpha)\mu_{g}$ (5)

The subscripts $l$ and $g$ denote gas and liquid phases respectively. The last
term in eq. (2) is the surface tension force, and $\sigma\kappa$ is the
pressure jump at the free surface. In our current algorithm, the jump
condition is enforced explicitly by representing the surface tension force as
a pressure gradient given by

$\displaystyle\sigma\kappa\mathbf{n}\delta\left(\mathbf{x}-\mathbf{x}_{f}\right)=-\nabla\tilde{p}$
(6)

where $\tilde{p}$ is the pressure solely due to the surface tension force at
the interface. An elliptic equation for $\tilde{p}$ is derived from continuity
and momentum equations,

$\displaystyle\nabla\cdot\left(\frac{\nabla\tilde{p}}{\rho}\right)=0$ (7)

The jump condition at the interface is enforced by using eq. (6) as a source
term in the solution of eq. (7). This ensures the exact difference in pressure
at the interface due to the surface tension. In our work, this elliptic
equation for $\tilde{p}$ is efficiently solved using a multigrid accelerated
red black SOR relaxation scheme implemented on a GPU. Using this method the
spurious velocities are seen to reduce to machine zero for a static bubble
with exact analytical curvature and to very small values when the curvature is
numerically computed. Admittedly, there are several methods to include the
surface tension force in the Navier-stokes (Renardy and Renardy [26], Sussman
et al. [27, 28], Gueyffier et al. [29]), however, in our experience the
current method more accurate for reduction of spurious velocities and
relatively easier to implement on GPUs.

###  2.2 Solution procedure

The following steps outline solution of the governing equations,

  1. 1.

Initialize the solution with initial liquid fraction ($\alpha^{n}$) and
velocities ($\mathbf{u}^{n}$)

  2. 2.

Calculate the density and viscosity at $n^{th}$ time step using,

$\displaystyle\rho^{n}=\alpha^{n}\rho_{l}+(1-\alpha^{n})\rho_{g}$
$\displaystyle\mu^{n}=\alpha^{n}\mu_{l}+(1-\alpha^{n})\mu_{g}$

  3. 3.

Solve the interface tracking equation (using $\mathbf{u}^{n}$) to obtain
$\alpha^{n+1}$

  4. 4.

Calculate the density ($\rho^{n+1}$) and viscosity ($\mu^{n+1}$) using
$\alpha^{n+1}$

  5. 5.

Calculate the interface curvature using the height function method

  6. 6.

Solve the pressure $(\tilde{p})$ using sharp surface force method

$\displaystyle\nabla\cdot\left(\frac{\nabla\tilde{p}}{\rho}\right)^{n+1}\hskip
14.22636pt\text{with}\leavevmode\nobreak\
\nabla\tilde{p}=-\left(\sigma\kappa\mathbf{n}\delta\left(\mathbf{x}-\mathbf{x}_{f}\right)\right)^{n+1}$

  7. 7.

Compute the pressure gradient term
$\left(\frac{\nabla\tilde{p}}{\rho}\right)^{n+1}$ at cell faces

  8. 8.

Compute the convection term $(\nabla\cdot\rho\textbf{uu})$ using geometry
construction method

  9. 9.

Compute the intermediate velocities at cell centers as,

$\displaystyle\rho_{c}^{n+1}\hat{\mathbf{u}}_{c}=\rho_{c}^{n}\mathbf{u}_{c}^{n}$
$\displaystyle-\Delta
t\sum_{f}\mathbf{u}_{c}^{n}\left(\sum_{k=1}^{2}\rho_{k}^{n}\alpha_{k}^{n}\right)A_{f}$
$\displaystyle+\Delta
t\sum_{f}\mu_{f}^{n}\left(\nabla\mathbf{u}^{n}+\nabla^{T}\mathbf{u}^{n}\right)_{f}\cdot\mathbf{n}_{f}+\rho_{c}^{n+1}\left(\frac{\nabla\tilde{p}}{\rho}\right)_{f\rightarrow
c}^{n+1}$

  10. 10.

Compute the intermediate velocities at cell faces as,

$\displaystyle\hat{\mathbf{u}}_{f}=\hat{\mathbf{u}}_{c\rightarrow f}$

  11. 11.

Compute pressure $(p)$ using,

$\displaystyle\nabla\cdot\left(\frac{\nabla
p}{\rho}\right)^{n+1}=\frac{\nabla\cdot\hat{\mathbf{u}}_{f}}{\Delta t}$

  12. 12.

Correct the velocity field for divergence free condition using,

$\displaystyle\mathbf{u}_{c}^{n+1}=\hat{\mathbf{u}}_{c}-\Delta
t\left(\frac{\nabla p}{\rho}+\frac{\nabla\tilde{p}}{\rho}\right)_{f\rightarrow
c}^{n+1}$

Where, c→f represents the linear averaging procedure to compute face value
(eg. $\rho_{f}$, $\mu_{f}$, $\mathbf{u}_{f}$) from its cell centered values.
The intricate details of the numerical algorithm is described in Kumar and
Vanka [17]. The entire algorithm has been programmed on the CUDA-Fortran
platform supported by the Portland Group (PGI) Fortran compiler and its detail
can be found in Kumar et al. [18].

##  3 GPU Implementation

This chapter gives a brief introduction to the Graphics Processing Unit (GPU),
describes its architecture and inner workings. The chapter also provides brief
details of the multiple GPU implementation and speed comparison.

###  3.1 Graphics Processing Units

The graphics processing unit (GPU) is a massively parallel device that can be
used for graphics generation and computing [30]. It is ubiquitously present in
every personal computers, gaming consoles and mobile phones. Until the last
decade, the GPUs have been primarily used for rendering the video games and
post processing high definition movies. The real-time rendering is a highly
parallel operation which requires multi-threaded high efficiency power horse
devices to perform the visualization operations. Along with the high-speed
rendering and video editing, they are also able to do fast single precision
and double precision arithmetics. With the introduction to Compute Unified
Device Architecture (CUDA) [31, 32], a parallel computing platform and
programming model the GPUs are also used in the scientific computing. Many new
generation supercomputers and data centers have dedicated GPU clusters. Using
CUDA, the computer programs can be written into ANSI C, ANSI C++ and Fortran.

Figure 1: A Pascal P100 GPU made by NVIDIA

Figure 1 shows a NVIDIA Pascal P100 GPU, which has a global memory of 16 GB
and a memory bandwidth of 720 GB/s. The GPU board contains over 15 Billion
transistors. It has a theoretical peak computation power of 9.3 TeraFLOPS for
single precision and 4.7 TeraFLOPS for double precision computations.

###  3.2 Multiple GPU Implementation

The computational performance of the GPU and adaptability of the CUDA
programming environment have led to widespread acceptance of the GPUs in the
scientific community. However, the exploration of fundamental sciences or
large scale engineering problems require significantly larger computational
resources than a single GPU can accommodate. It is true that the GPU hardware
has grown significantly over the time a relatively smaller period of time. For
example, the global memory per GPU has increased from 4 GB in 2009 to 24 GB in
2015. The number of CUDA core per GPU have grown from 192 in 2009 to 3584 in
2016. But it is still not sufficient to study a conduct direct numerical
simulation of bubbly flow in a square duct. The obvious solution to this
problem is to utilize multiple GPUs simultaneously. The multiple GPU (multi-
GPU) implementations heavily borrows the programming experiences from the
Massage Passing Interface (MPI) and the OpenMPI programming environments.

There are three main sections of the multi-GPU implementations:

  1. 1.

Domain decomposition

  2. 2.

Communication between adjacent domains

  3. 3.

Optimization of communication and computation times

####  3.2.1 Domain Decomposition

As name suggests, the domain decomposition refers to splitting of a larger
computational domain into smaller subdomains which can fit on a single GPU. In
case of three-dimensional Cartesian domain, the decomposition can be performed
in several ways e.g. along $x$, $y$ or $z$ directions or $x$ and $y$
directions or $y$ and $z$ directions or $x$, $y$ and $z$ directions. Any of
these direction can be chosen, however, a load balancing must be performed for
its optimization. In our current implementation, we have chosen to split the
domain along $z$ direction. This minimizes the communication time and reduces
the programming difficulty. Figure 2 shows an illustration of one-directional
domain decomposition along the horizontal direction of a two-dimensional
computational domain. The yellow and gray colored cells represent the interior
and boundary cells respectively. It can be noticed from figure that equal
number of computational cells has been prescribed in each subdomain. This is
enforced to ensure optimum load balancing for each GPU.

Figure 2: Illustration of Domain decomposition for 4 GPUs

The domain decomposition alleviates the shortcomings of smaller global memory.
But, it creates virtual boundaries between two neighboring domains. Hence, a
mechanism to pass the information between adjacent domain/GPU is required. To
accommodate the message passing between neighbors one layer is added on each
side of the domain and these cells are referred as halo cells. Figure 3 shows
the subdomains on each GPUs and their respective data structures. The LO and
RO are the left and right outer cells and used to facilitate the data transfer
between neighboring domains. For example, the right outer cells of GPU0
mirrors the left interior of GPU1 and left outer cells of GPU1 mirrors the
right interior of the GPU0. The LI and RI are the left and right interior
cells which is useful to facilitate physical boundary conditions on GPU0 and
GPU3 respectively. The left and right interior cells are also used to transfer
the data between neighboring cells. Although The left and right outer cells
are not needed on GPU0 and GPU3 respectively, they are added to maintain
uniform data structures across the GPUs. It should be also noted that the
cells numbering on subdomains are different from that on the global domain.

Figure 3: Additional layers at the domain boundaries

####  3.2.2 Communication between Neighboring Domains

Figure 4 illustrates the data transfer procedure between the two nodes
connected via network. It is a three step process. In the first step, the data
is copied from the GPU0 memory to the CPU0 memory on Node 0 via PCIe
interface. In the second step, the data is sent from the Node 0 to Node 1 via
InfiniBand [33]. The Message Passing Interface (MPI) [34, 35] is used to
transfer the data over the network. In the last step, the data is copied from
the CPU1 memory to the GPU1 memory on Node 1. In these three steps, the data
is copies from the GPU 0 to the GPU 1. This process is performed co-currently
on every Node to minimize the down/waiting time. Once the transfer from left
to right ($0\rightarrow 1$, $1\rightarrow 2$, $2\rightarrow 3$ etc.) has taken
place, the transfer from right to left ($0\leftarrow 1$, $1\leftarrow 2$,
$2\leftarrow 3$ etc.) is initiated. The total time taken to complete the data
copying process is referred as communication time.

Figure 4: Data transfer process between two nodes

On a shared memory machine where a multi-core CPU and multiple GPUs are
present, the data copying process is less cumbersome. During the course of
calculation, each GPU is assigned to a particular CPU core. The copying
process is performed using same routine as explained earlier. The first and
last steps of the data transfer is same to that of the distributed systems. In
the second step, the data is passed from one CPU to another via pointers.
Although absence of data copying over network reduces the overall
communication time the reduction is relatively less prominent. This technique
to balance computation time over communication time has also been applied in
our Lattice Boltzmann implementation [36, 37, 38]

####  3.2.3 Multi-GPU Speedup

The performance of multi-GPU implementation has been tested on a four GPU
shared memory workstation as well as on Blue Waters (BW) supercomputer which
is an distributed memory cluster. It should be noted that processing speed of
in-house workstation CPU is better than that of Blue Water CPU. A multi-CPU
version code has also been developed for CPU vs GPU speedup comparisons. A 3D
lid-driven cavity problem is selected to test the performance of CUFLOW. The
governing equations are solved on a grid of $128\times 128\times 512$ or
approximately 8 million control volumes.

Table 1: Specification of the hardwares used |  In-house workstation |  Blue waters
---|---|---
# of nodes |  1 |  4224
Node CPU |  Xeon E5-2650v2 |  AMD 6276
|  2.60 GHz, 8 cores |  2.3 GHz, 16 cores
GPU/Node |  4 $\times$ NVIDIA Tesla C2075 |  1 $\times$ NVIDIA Tesla K20x
|  5 GB each, 575 MHz |  6 GB each, 732 MHz
Figure 5: Time taken by multi-CPU and multi-GPU CUFLOW

Figure 5 shows the time taken to complete one time step of the simulation. It
can noticed that the CPU version of the CUFLOW takes longer time on the Blue
Water supercomputer than on the in-house workstation. However, the GPU version
takes lesser time on Blue water. The multi-GPU code on in-house workstation
reduces the computation time by a factor of 3.2 between one GPU and four GPUs
whereas on the Blue Water supercomputer the factor is approximately 2. This is
primarily due to increased communication time overhead.

##  4 Computational details

The simulations presented here are performed for an upward bubbly flow in a
square duct. In the $x\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ y$
directions, the flow is wall-bounded where no-slip and no-penetration
conditions are applied. In the $z$ or streamwise direction the flow is
considered periodic. The gravitational acceleration acts in the downward or
negative $z$ direction, and the contribution due to hydrostatic pressure on
the continuous phase is removed with the addition of $(\rho_{0}\mathbf{g})$,
where $\rho_{0}=(1-\bar{\phi})\rho_{l}+\bar{\phi}\rho_{l}$
$(\bar{\phi}:\text{average void fraction})$ to Navier-stokes equations. This
was employed in order to prevent a uniform downward acceleration of the
continuous phase. A constant mean pressure gradient based on bulk Reynolds
number $Re_{b}=\rho_{l}w_{b}W\mu_{l}^{-1}$ of unladen flow has been applied
along positive $z-$direction. It has been kept constant between unladen and
laden flows, hence the bulk velocity of laden flow is expected to increase as
the average density within the domain will decrease. The bulk Reynolds number
of unladen flow simulation is set to 5300, so that sustainable turbulence is
realized. The corresponding frictional Reynolds number
$Re_{\tau}=\rho_{l}u_{\tau}W\mu_{l}^{-1}$ is approximately 360.

Two simulations of turbulent flow: (a) unladen flow, and (b) bubble laden flow
are performed in the present study. Both simulations are performed on a
Cartesian grid of $192\times 192\times 768$ control volumes with a uniform
grid spacing in $x$, $y$, and $z$ directions. The extent of a square duct in
$x$, $y$, and $z$ directions is $48mm\times 48mm\times 192mm$. We note that
the canonical height $(H)$ used to understand unladen turbulence flow in ducts
and channels is $H=2\pi W$ [39, 40] or higher, however, for bubble laden flow,
it has been observed that the correlation length of two-phase flow is smaller
than that of single phase flow [8, 41, 10, 42]. Hence, for computational
efficiency we have selected a smaller duct length (four times the duct width).
Admittedly, turbulence quantities of unladen flow are not expected to match
exactly with results available in literature [40]. However, since the primary
focus of this work is on understanding the mechanisms for capturing and
redistribution of the bubble by turbulence and modification of the velocity
structures, root mean square velocities and Reynolds stresses by the dispersed
phase this domain size is adequate.

The unladen flow is started with a perturbed solenoidal velocity and
integrated till a stationary flow field is realized. The initial bulk velocity
and fluctuations are given by,

$\displaystyle
w(x,y,z)=\begin{cases}w_{b}\left(\frac{2y}{W}\right)^{1/7}&y/W\in[0,0.5]\\\
w_{b}\left(2-\frac{2y}{W}\right)^{1/7}&y/W\in[1,0.5]\end{cases}$ (8)
$\displaystyle u^{\prime}(x,y,z)=0.08w_{b}\cdot\text{cos}\left(\frac{n\pi
x}{W}\right)\frac{W}{H}\left[\text{cos}\left(\frac{n\pi
z}{H}\right)\text{sin}\left(\frac{n\pi
y}{W}\right)+(\text{sin}\left(\frac{n\pi
z}{H}\right)\text{cos}\left(\frac{n\pi y}{W}\right)\right]$ (9) $\displaystyle
v^{\prime}(x,y,z)=0.08w_{b}\cdot\text{sin}\left(\frac{n\pi
z}{H}\right)\text{sin}\left(\frac{n\pi y}{W}\right)\text{sin}\left(\frac{n\pi
x}{W}\right)$ $\displaystyle
w^{\prime}(x,y,z)=0.08w_{b}\cdot\text{sin}\left(\frac{n\pi
z}{H}\right)\text{sin}\left(\frac{n\pi y}{W}\right)\text{sin}\left(\frac{n\pi
x}{W}\right)$

The bulk velocity $w_{b}$ of unladen flow is calculated from the bulk Reynolds
number as $w_{b}=Re_{b}\mu_{l}\rho_{l}^{-1}W^{-1}$. The one-seventh law of
turbulent channel flow with $8\%$ divergence free velocity fluctuations is
prescribed as the initial condition for velocity field. Then the governing
equations are integrated with the initial conditions and after a sufficient
time, $t=60T_{b}$, where $T_{b}=H/w_{b}$ is the bulk time, the effects of
initial condition are seen to diminish and a stationary turbulence field is
realized. At this point, we initiate the calculation of mean variables and
collect them until fully developed turbulence is realized and a stationary
mean velocity profile is achieved. The fully developed turbulence (or
stationary mean flow) is achieved by approximately $t=145T_{b}$. Turbulence
statistics such as rms velocities $(\sqrt{\left\langle
u^{\prime}u^{\prime}\right\rangle},\sqrt{\left\langle
v^{\prime}v^{\prime}\right\rangle},\sqrt{\left\langle
w^{\prime}w^{\prime}\right\rangle})$ and Reynolds stresses $(\left\langle
u^{\prime}u^{\prime}\right\rangle,\left\langle
v^{\prime}v^{\prime}\right\rangle,\left\langle
w^{\prime}w^{\prime}\right\rangle,\left\langle
u^{\prime}v^{\prime}\right\rangle,\left\langle
u^{\prime}w^{\prime}\right\rangle)$ are collected for an additional 55 bulk
times (or until $t=200T_{b}$). The results of single phase calculations are
presented in section 5.1.

The turbulent flow field from the solution of unladen flow is taken as the
initial condition for laden flow calculation. At $t=200T_{b}$, 864 mono-
dispersed bubbles are uniformly placed in the computational domain. Initial
velocities inside the bubbles are set to that of the liquid velocities at
respective spatial locations. Based on the total volume of bubbles, the void
fraction of the laden flow is approximately $\bar{\phi}=0.82\%$. Since void
fraction is less than $2\%$, the mixture can be considered dilute and
hydrodynamic interactions among bubbles are expected to be negligible. The
bubbles are arranged in a regular array of $6\times 6\times 24$ bubbles in
$x,y,\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ z$ directions
respectively. Although they are placed in a regular array, their initial
configuration is expected to have no effect on the stationary-state results.
Bunner and Tryggvason [43, 44] have explored the effects of the initial
configuration of bubbles and found negligible effect on the final solution.

Figure 6 shows the illustration of computational domain and the initial
configuration of the top and bottom layers of 36 bubbles (the other 22 layers
are not shown in the picture, but they are present in the simulations). The
liquid and gas densities are 1000 and 1.2 kg/m3 respectively, which results in
a density ratio of $\rho_{l}/\rho_{g}=833$. The dynamic viscosity of the
liquid and gas are $0.001$ and $1.8\times 10^{-5}$ Pa.s, respectively. The
surface tension $(\sigma)$ between liquid and gas is $0.072$ N/m. The diameter
of each bubble is $2$ mm. The Bond number $Bo=\rho_{l}gd^{2}/\sigma$ (or
Eötvos number $Eo=\Delta\rho gd^{2}/\sigma$) and Morton number
$Mo=g\mu_{l}^{4}/\rho_{l}\sigma^{3}$ are approximately 0.55 and $2.6\times
10^{-11}$, respectively. Based on the combination of Bond and Morton numbers,
the expected terminal shape of a bubble rising in a quiescent liquid column is
nearly spherical [45], and it is expected to follow a vertical straight path
during its rise. The choice of smaller bubble is to explore the effects of
spherical bubbles. However, due to local shear and turbulence, a small degree
of deformation as well as a non-rectilinear rise path is expected.

Figure 6: Computational domain and initial bubble locations

##  5 Result and discussion

###  5.1 Validation: Unladen Duct Flow

We present results of a validation study for predictions of single phase
turbulence and compare them with results in literature. The single phase
turbulence predictions are verified qualitatively and quantitatively with the
works of Gavrilakis [46], Huser and Biringen [47], Broglia et al. [48], Zhang
et al. [49]. The frictional Reynolds number $Re_{\tau}=360$ of the current
work is intermediate to Gavrilakis [46], Huser and Biringen [47], Zhang et al.
[49], and it is same as Broglia et al. [48], a large eddy simulation, hence,
only qualitative comparisons can be made for some turbulence quantities. The
algorithm used for predicting turbulent flows in square ducts has also been
validated by Winkler et al. [50, 51], Chaudhary et al. [52, 53, 54], Shinn
[55]. The presented quantities are averaged over four statistically invariant
transformations: time, $z$-direction, to account for periodicity, and two
central planes: $x$ and $y$, to enforce quadrant symmetry.

Figure 7: Averaged streamwise velocity $(w)$ y-wall bisector

Figure 7 shows the mean streamwise velocity distribution along the wall
bisector. It has been scaled by the centerline streamwise velocity $(w_{c})$
which is $w_{c}=0.1464$. It can be noticed that the current prediction of
streamwise velocity profile is very close to the predictions made by
Gavrilakis [46] for $x>0.0072\leavevmode\nobreak\ mm$ (or $x>0.3W)$. For
$x<0.3W$, the predicted mean velocity is approximately 1-2% higher. It should
be noted that they used a stretched grid of $127\times 127$ nodes in the
cross-section and their first grid node was at $\Delta x^{+}=0.45$, whereas
the current results are obtained using a uniform grid of $192\times 192$ nodes
with the first node at $\Delta x^{+}=1.875$. Hence, a small difference in the
near-wall statistics is expected between these two studies. A similar
qualitative agreement can also be observed between current and Broglia et al.
[48] predictions. The current results are nearly identical for $x\leq 0.1W$
and are within 1-2% for $x>0.1W$. It is also noteworthy that point-wise data
for comparison purposes were extracted from the original articles using Matlab
imaging tool and manual hand picking, hence a small human error can also be
expected.

Figure 8 shows the streamwise $(w_{rms})$ and spanwise $(v_{rms})$ velocity
fluctuations along the y-wall bisector for $x^{+}\leq 100$. A qualitative
agreement between current predictions and Zhang et al. [49] prediction at
$Re_{\tau}=300$ and Huser and Biringen [47] at $Re_{\tau}=600$ can be noticed
for both $(w_{rms})$ and $(v_{rms})$ from the figs. 8(a) and 8(b),
respectively. In direct comparison with $Re_{\tau}=360$, it is evident that
$v_{rms}$ from our prediction is in close agreement with that of Broglia et
al. [48], whereas the streamwise rms velocity, $w_{rms}$, is identical to
current work for $x^{+}\leq 18$ and higher for $x^{+}>18$. The higher value of
$w_{rms}$ in Broglia et al. [48] is the artifact of a larger streamwise grid
spacing. A similar effect of grid resolution has also been observed by Zhang
et al. [49].

(a)

(b)

Figure 8: Root mean square velocity fluctuations along the wall bisector: (a)
$w_{rms}$ and (b) $v_{rms}$

Another way to measure the solution’s accuracy is to compare the ratio of
centerline streamwise and bulk velocities, $(w_{c}/w_{b})$, and the ratio of
frictional and bulk velocities, $(w_{\tau}/w_{b})$ to those reported in the
literature. The values of $(w_{c}/w_{b})$ and $(w_{\tau}/w_{b})$ are 1.33 and
0.678, respectively which are within 0.3% of that reported in the literature
[40, 48]. These simulations validate the grid used and provide data for
comparisons with bubble laden flow.

###  5.2 Bubble Laden Flow

Now, we present the results of study to analyze the effects of bubbles on
liquid phase turbulence and discuss interaction between bubbles and
turbulence. We discuss the findings of unladen flow and compare them with the
bubble laden flow results. The instantaneous flow discussions are presented
first, followed by a discussion on the turbulence statistics. Most of the
turbulence quantities are calculated and discussed near the wall region, along
the corner bisector and along the wall bisector. These locations are shown in
fig. 9.

Figure 9: Locations along which the quantities are presented

####  5.2.1 Bubble Distribution

Figure 10 shows the initial placement of bubbles in the unladen flow. The
isometric, front and top views are shown in Figures 10(a), 10(b) and 10(c).
Furthermore, in order to maintain a better aspect ratio of the figures, only
half of the computational domain is shown in the figs. 10(a) and 11(a).

(a) Isometric view

(b) Front view

(c) Top view

Figure 10: Initial configuration of bubbles in the periodic domain

Figure 11 shows the transient distribution of bubbles at $t_{laden}=40T_{b}$
where $t_{laden}$ is the total time for which the bubble laden flow
calculation has been conducted. The isometric view given in fig. 11(a) shows
that the bubble are randomly distributed in the domain. However, a clearer
idea of bubble distribution is gathered from the front and top views shown in
figs. 11(b) and 11(c) respectively.

Comparing figs. 10 and 11, we can state that the bubbles are significantly
affected by the flow as they have moved away from their initial positions.
From the top view it can be noticed that a majority of bubbles have migrated
towards the wall and formed a layer. These layers of bubbles are very close
(less than a control volume) to the wall111We have used a complete wetting
boundary condition for the sidewalls which ensure a thin layer of liquid at
walls. Due to this the bubbles are able to get very close the wall without
making a contact with them. as seen in both the front and top views. Away from
the first layer near the wall, a cluster of bubbles are observed, and they are
located in the regions between the wall and the half distance from the duct
center. Furthermore, the cluster is sparse and skewed i.e. a larger number of
bubbles are present near the corner region than that near the wall bisector
region. The presence of a bubble layer near the wall and a cluster of bubbles
away from the wall are clear indications that the current scenario is wall
wall-peaking with a secondary intermediate-peak.

(a) Isometric view

(b) Front view

(c) Top view

Figure 11: Instantaneous configuration of bubbles as seen in (a) isometric,
(b) front and (c) top views at $t=274$ s or after 40 bulk time units after
introduction of bubbles in the unladen flow

The ensemble-averaged distributions of the void fraction
$\left\langle\phi\right\rangle$ along wall bisector, wall and corner bisector
are shown in figs. 12(a), 12(b) and 13(a). They are symmetric about duct
center and wall bisector therefore only half of the domain are shown in the
figures. It can be noticed from figs. 12(a) and 13(a) that there are two peaks
along the wall and corner bisectors. The observation is in the direction as
that observed in fig. 11. The thickness of the peak at these locations are
approximately 2 mm, or one bubble diameter. Furthermore, the peak value along
the corner bisector is higher that along the wall; therefore, the bubble
density near the corner is larger than that near the wall bisector.

In the wall direction there are multiple peaks and valleys which indicate that
there are a large number of bubbles close to the wall, and there is a
significant gap between these bubbles. The height and width of the first peak
is similar to that along the corner bisector. Rest of peaks are smaller and
wider than the first peak, hence, the bubble density along the wall is highest
in the corner region, and it is nearly constant away from the corner region.

(a) Along wall bisector

(b) Along wall

(a) Along corner bisector Figure 13: Ensemble-averaged void fraction profiles
along (a) wall bisector, (b) wall and (c) corner bisector at $t=274$ s

####  5.2.2 Instantaneous Flow Structures

The three-dimensional fluctuation of the instantaneous streamwise velocity
$(w^{\prime}_{3d})$ can give a qualitative picture of the flow field. The
expression for $w^{\prime}_{3d}$ is given by

$w^{\prime}_{3d}(x,y,z,t)=w(x,y,z,t)-\left\langle w\right\rangle_{z}(x,y,t)$
(10)

The instantaneous streamwise velocity, $w(x,y,z,t)$, is averaged along the
periodicity direction to get $\left\langle w\right\rangle_{z}(x,y,t)$. The
average is performed only in the regions where liquid is present. The iso-
surfaces and contours of instantaneous streamwise velocity fluctuations
$w^{\prime}_{3d}$ are shown in footnote 2.

Figure 14(a) shows the instantaneous streamwise flow structures in the unladen
flow. It can be noticed that the instantaneous flow is dominated by elongated
structures in the streamwise direction, and they are mostly located in the
near wall regions and extend nearly half of duct in the streamwise direction.
The magnitude of instantaneous fluctuations $w^{\prime}_{3d}$ is approximately
half of bulk velocity.

In the case of laden flow, the instantaneous flow consists of both positive
and negative fluctuations, but the long elongated structures observed in
unladen flow split into a set of smaller and slimmer structures as shown in
fig. 14(b). The wall region is dominated by relatively longer instantaneous
structures, whereas the core region is dominated by smaller structures. The
smallest size of structures in the core region is comparable to the bubble
size. These structures, in both core and wall regions, are more oscillatory
than the ones observed in the unladen flow. The magnitude of fluctuations are
approximately ten times larger than of unladen flow.

(a)

(b)

Figure 14: Instantaneous three-dimensional velocity fluctuations for (a)
unladen flow and (b) bubble laden flow. Displayed are the iso-surface222It
should be noted that the centerline velocity of unladen flow is approximately
2.5 times smaller than that of bubble laden flow. Hence, the value of iso-
surfaces in figs. 14(a) and 14(b) are not same.of
$w^{\prime}_{3d}/w_{c}=0.15$, the vertical distance between planes used for
contour plots is $W/2$.

####  5.2.3 Instantaneous Streamwise Velocity

Instantaneous flow field of unladen flow in the duct cross-section is shown in
fig. 15(a). The velocities are averaged along the streamwise direction to
increase the sample size. The figure shows the contour of instantaneous
streamwise velocity overlaid with instantaneous in-plane velocity vectors.
Significant fluctuations in the velocity can be observed from the contours of
the velocity, which is indicative of the turbulent flow. The presence of the
in-plane velocity vectors are also a sign of turbulent flow. It can be noticed
that there is a large variation in the magnitude of the in-plane velocities,
indicated by the length of vectors, at similar locations in different
quadrants. The phenomenon is due to the dynamic nature of the turbulent flow.

Figure 15(b) shows the instantaneous streamwise velocity of the bubble laden
flow. It can be noticed that the streamwise velocity fluctuations are much
larger compared to the unladen flow, and they are as high as bulk velocity.
Furthermore, the maximum fluctuations are observed in the corner and wall
regions where bubbles are located. The streamwise velocity fluctuations in the
vicinity of the bubbles are comparable to the terminal velocity of the bubble.
The bulk velocity of the bubble laden flow is close to 0.322 m/s which is
nearly three times higher than that of unladen flow. The increase of bulk
velocity is a combined result of the upward motion of bubbles (and subsequent
dragging of the surrounding liquid) and significant reduction in viscosity
near the walls due to the preferential accumulation of bubbles. The cross-
stream velocities are also higher than their unladen flow counterparts, and
they are located near wall and corner regions.

(a)

(b)

Figure 15: Contours of instantaneous streamwise velocity overlaid with cross-
stream velocity vectors for (a) unladen flow and (b) bubble laden flow

####  5.2.4 Ensemble-averaged Distribution of Streamwise and Spanwise Velocity

The contours of ensemble-averaged333The ensemble average of r.m.s and
turbulence quantities ($\left\langle u\right\rangle$,$\left\langle
v\right\rangle$,$\left\langle w\right\rangle$, $\left\langle
u^{\prime}u^{\prime}\right\rangle$, $\left\langle
v^{\prime}v^{\prime}\right\rangle$, $\left\langle
w^{\prime}w^{\prime}\right\rangle$, $\left\langle
u^{\prime}v^{\prime}\right\rangle$, $\left\langle
u^{\prime}w^{\prime}\right\rangle$, $\left\langle\phi\right\rangle$) are
performed over time, $z$-direction, and $x$ and $y$ central planes. Then the
averaged values are reflected back to other three quadrants and the
appropriate sign of the quantities are maintained. Hence, the quantities are
expected to be quadrant symmetric. streamwise velocity of unladen and laden
flows in the duct cross-section are given in figs. 16(a) and 16(b). The
contours are overlaid with the ensemble-averaged in-plane velocity
($\left\langle u\right\rangle$,$\left\langle v\right\rangle$) vectors. It can
be noticed that there are two counter-rotating vortices in each quadrant and a
total of eight vortices in the duct cross-section. They are the results of the
mean secondary flow generated by anisotropy of turbulent stresses. They are
responsible for the transfer of high energy eddies from the duct center to the
walls. The mean secondary flow transports the fluid momentum from the bulk
region to the corners along their bisectors, and then back to the bulk region
along wall bisectors. The magnitude of mean secondary velocities reach a
maximum along the diagonal of the duct.

Comparing unladen and bubble laden flows, we notice that the streamwise
velocity increases with the introduction of bubbles. The centerline velocity,
$w_{c}$, increases from $0.146$ m/s for unladen flow to $0.325$ m/s for bubble
laden flow. The velocity increases everywhere in the cross-section and not
just at the centerline. It can also be observed that the streamwise velocity
of bubble laden flow is flatter than that of unladen flow and the velocity
contour lines in the corner region are rounder in laden flow case than that in
unladen flow.

(a)

(b)

Figure 16: Contours of the mean streamwise velocity overlaid with mean cross-
stream velocity vectors for (a) unladen flow and (b) bubble laden flow

Figures 17(a) and 17(b) show the mean in-plane streamlines in the cross-
sectional plane. It can be noticed that the introduction of bubbles in unladen
flow does not change direction of the secondary flows, i.e. fluid from the
bulk region to corners are transferred along the corner bisectors and then
back to the bulk region along wall bisectors. However, the introduction of
bubbles makes the flow a bit more complex and introduces additional counter-
rotating vortices near the corner bisectors. They are absent in the unladen
flow, and streamlines reach the wall before they bend. The additional vortices
in the laden flow affect the streamline velocity $w$ distribution in the
corner region, as it is reflected in the form of a pointed counter line in
fig. 16(b). In the unladen flow, counter-rotating vortices are nearly
symmetric about the corner bisectors, whereas in the laden flow, they are
clearly asymmetric. It can also be observed that the secondary flow
streamlines in the absence of bubble are smooth and free of oscillations,
whereas they are tortuous and twisted when bubbles are introduced. This may be
because the bubble distribution is not fully symmetric and/or the flow can
never reach full symmetry. It is also possible that the flow has not reached a
stationary state and more sampling of the results are required.

(a)

(b)

Figure 17: Contours of mean streamwise velocity overlaid with mean cross-
stream stream functions for (a) unladen flow and (b) bubble laden flow

(a)

(b)

Figure 18: Contours of mean $x$-direction velocity in the lower left cross-
section of the duct: (a) unladen flow and (b) bubble laden flow

Contours of $x$-direction spanwise velocity of unladen and laden flows are
shown in figs. 18(a) and 18(b). The maximum magnitude of $\langle u\rangle$ is
approximately 2% of streamwise bulk velocity $(w_{b})$ for both unladen and
laden flows. The unladen flow prediction is in accordance with the data
available in literature [40, 48]. It can be noticed that in the corner region,
$\langle u\rangle$ in laden flow case is approximately four times higher than
that in unladen flow, and the highest value is observed closer to the corner
in laden flow, whereas in unladen flow the highest value is observed at a
farther distance from the corner. It can also be seen that $x$-velocity is
nearly similar in both laden and unladen flows in the regions away from the
corner. This higher value of $\langle u\rangle$ near the corner region is
associated with the presence of bubbles.

(a)

(b)

Figure 19: Contours of mean $y$-direction velocity in the lower left cross-
section of the duct: (a) unladen flow and (b) bubble laden flow

Contours of $y$-direction spanwise velocity are shown in figs. 19(a) and
19(b). In the case of unladen flow, the maximum magnitude of $\langle
v\rangle$ is close to 2% of the streamwise bulk velocity444The bulk velocity
of the unladen flow is 0.110 m/s $(w_{b})$. This observation is similar to
that for $\langle u\rangle$. In the case of laden flow, the maximum magnitude
of $\langle v\rangle$ is smaller than that of $\langle u\rangle$. This is
somewhat surprising as we would expect both $\langle u\rangle$ and $\langle
v\rangle$ to be similar as there is no source of directional preference and
can be attributed to a lower sampling time of the simulation. The highest
value of $\langle v\rangle$ in laden flow is three times higher than that in
the unladen flow. Both location of the highest value and distribution of
$\langle v\rangle$ in the bulk region is similar in unladen and laden flows.

It can also be noticed from figs. 18(a) and 19(a) that $\langle u\rangle$ is
an approximate reflection of $\langle v\rangle$ about corresponding diagonal.
The slight asymmetry between $\langle u\rangle$ and $\langle v\rangle$ is a
mere indication of the quadrant, not octagonal, averaging. It should also be
noted that the range of $x$ and $y$ velocities are nearly identical and in
accordance with physics of the problem. However, in the case of laden flow, no
such parallels can be drawn between figs. 18(b) and 19(b).

Although the mean secondary flow velocities are two orders of magnitude
smaller than streamwise velocity, it has a significant impact on the
distribution of streamwise velocity. As, it can be observed from fig. 16(a)
that the $w$-contour lines near the corner are pointed and bulged away from
the duct corner. The mechanism is related to the two counter-rotating vortices
present near the corner.

####  5.2.5 Reynolds stresses

The Reynolds stresses are presented along all three locations, namely wall
bisector, wall and corner bisector as mentioned in fig. 9. The data points are
extracted using a visualization software; thus we expect a small interpolation
errors in the results. Although the laden flow has not reached a stationary-
state, we can analyze the trends of the Reynolds stresses and compare them
with unladen flow to understand the effects of bubbles.

####  5.2.6 Along Wall Bisector

The normal components $\langle u^{\prime}u^{\prime}\rangle$, $\langle
v^{\prime}v^{\prime}\rangle$, and $\langle w^{\prime}w^{\prime}\rangle$ of
Reynolds stress along the wall bisector for both unladen and bubble laden
flows are given in fig. 20. All the components of Reynolds stress are scaled
with the square of unladen flow bulk velocity. It can be noticed from the
figs. 20(a), 20(b), 20(c), 20(d), 20(e) and 20(f) that $\langle
u^{\prime}u^{\prime}\rangle/w_{b}^{2}$, $\langle
v^{\prime}v^{\prime}\rangle/w_{b}^{2}$, $\langle
w^{\prime}w^{\prime}\rangle/w_{b}^{2}$ are zero at the wall and symmetric
about $x=0.024$ or $W/2$ which is expected as the quadrant symmetry has been
enforced.

In the case of unladen flow, all three normal components initially increase
with distance from the wall to reach a maximum and then decrease to reach a
minimum. The location of the maximum is closer to the wall, and the minimum is
at center. It can also be noticed that the rate of increase of these
quantities is higher than the rate of their decrease. The corresponding
maximums of $\langle u^{\prime}u^{\prime}\rangle/w_{b}^{2}$, $\langle
v^{\prime}v^{\prime}\rangle/w_{b}^{2}$, $\langle
w^{\prime}w^{\prime}\rangle/w_{b}^{2}$ are observed at $x=0.0041$, 0.0054 and
0.0018 or $x=0.085W$, $0.113W$ and $0.037W$, respectively. The maximum values
of $\langle w^{\prime}w^{\prime}\rangle/w_{b}^{2}$ and $\langle
v^{\prime}v^{\prime}\rangle/w_{b}^{2}$ are roughly 10 and 1.5 times higher
than that of $\langle u^{\prime}u^{\prime}\rangle/w_{b}^{2}$. It can also be
noticed that the ratio between maximum and centerline values of $\langle
u^{\prime}u^{\prime}\rangle$, $\langle v^{\prime}v^{\prime}\rangle$, and
$\langle w^{\prime}w^{\prime}\rangle$ are 2, 3 and 11 respectively. These
values have also been predicted earlier in Madabhushi and Vanka [40].

In the case of laden flow, there are multiple local maximum and minimum for
all three normal components. It can be seen from figs. 20(b), 20(d) and 20(f)
that the maximum is always observed near the wall, and the minimum of $\langle
u^{\prime}u^{\prime}\rangle$ and $\langle v^{\prime}v^{\prime}\rangle$ are
observed at the center. The minimum of $\langle w^{\prime}w^{\prime}\rangle$
is observed between the wall and duct center. It can also be noticed that the
values of all three quantities at all locations are considerably higher in the
presence of bubbles . For example, at the center $\langle
u^{\prime}u^{\prime}\rangle/w_{b}^{2}$ in unladen case is 0.0024, whereas in
laden flow, it is approximately 0.0065.

Similar to unladen flow, the maximum value of $\langle
w^{\prime}w^{\prime}\rangle/w_{b}^{2}$ and $\langle
v^{\prime}v^{\prime}\rangle/w_{b}^{2}$ is approximately 8 and 3 times higher
than $\langle u^{\prime}u^{\prime}\rangle/w_{b}^{2}$ for bubble laden flow as
well. Trend that the first peak of $\langle
w^{\prime}w^{\prime}\rangle/w_{b}^{2}$ is closer to the wall than that of
$\langle v^{\prime}v^{\prime}\rangle/w_{b}^{2}$ and the peak of $\langle
v^{\prime}v^{\prime}\rangle/w_{b}^{2}$ is closer to the wall than that of
$\langle u^{\prime}u^{\prime}\rangle/w_{b}^{2}$ is consistent between unladen
and bubble laden flows. However, there is one key difference in shape of the
$\langle w^{\prime}w^{\prime}\rangle/w_{b}^{2}$. In the case of unladen flow,
the curve is ‘M’ shaped, whereas in the case of bubble laden flow, the curve
is ‘W’ shaped.

(a) $\langle u^{\prime}u^{\prime}\rangle$ component for unladen flow

(b) $\langle u^{\prime}u^{\prime}\rangle$ component for bubble laden flow

(c) $\langle v^{\prime}v^{\prime}\rangle$ component for unladen flow

(d) $\langle v^{\prime}v^{\prime}\rangle$ component for bubble laden flow

(e) $\langle w^{\prime}w^{\prime}\rangle$ component for unladen flow

(f) $\langle w^{\prime}w^{\prime}\rangle$ component for bubble laden flow

Figure 20: Reynold stress components: (a-b) $\langle
u^{\prime}u^{\prime}\rangle$, (c-d) $\langle v^{\prime}v^{\prime}\rangle$ and
(e-f) $\langle w^{\prime}w^{\prime}\rangle$

Figures 21(a), 21(b), 21(c) and 21(d) show the $\langle
u^{\prime}v^{\prime}\rangle$ and $\langle u^{\prime}w^{\prime}\rangle$
components of Reynolds stress along the wall bisector for both unladen and
bubble laden flows. It can be noticed that these two quantities are asymmetric
about the center or $x=W/2=0.024$. For both unladen and laden flows, $\langle
u^{\prime}v^{\prime}\rangle$ is nearly zero along the wall bisectors. This
signifies the very small fluctuations of $x$ and $y$ velocities along the wall
bisector.

(a) $\langle u^{\prime}v^{\prime}\rangle$ component for unladen flow

(b) $\langle u^{\prime}v^{\prime}\rangle$ component for bubble laden flow

(c) $\langle u^{\prime}w^{\prime}\rangle$ component for unladen flow

(d) $\langle u^{\prime}w^{\prime}\rangle$ component for bubble laden flow

Figure 21: Reynold stress components: $\langle u^{\prime}v^{\prime}\rangle$
and $\langle u^{\prime}w^{\prime}\rangle$

It can be noticed from figs. 21(c) and 21(d) that the $\langle
u^{\prime}w^{\prime}\rangle$ initially decreases with the distance to reach a
minimum and then increases to reach a maximum at the same location on the
other side of the duct. In the unladen flow, the minimum and maximum are
observed at $x=0.0044$ or $0.091W$ and $x=0.909W$ respectively. The minimum
and maximum values are -0.004 and 0.004 respectively, which are approximately
same as $\langle u^{\prime}u^{\prime}\rangle/w_{b}^{2}$. In the case of laden
flow, the minimum and maximum are observed at $x=0.0044$ and $x=0.0435$
respectively which are also at the same locations where they were observed in
the unladen flow. The shape of $\langle u^{\prime}w^{\prime}\rangle$ in the
central region is also similar to each other as their approximate slopes are
0.2042 and 0.2837 respectively in unladen and laden flows.

From figures of unladen and laden flow Reynolds stresses, we can make two
qualitative conclusions: (1) Reynolds stresses increase when bubbles are
introduced in the flow and (2) their maximum modification are observed near
the wall. Due to the added effects of the buoyancy forces, bubbles travel with
higher speed than the surrounding liquid in the upward flow, hence it
introduces additional fluctuations in the flow; thus, we observe an increase
in velocity fluctuations everywhere in the domain. The second observation is
related to the migration of bubbles toward the wall region. As they
concentrate near the wall, they increase the velocity fluctuations.

##  6 Conclusions

In the present study, the dynamics of a bubble swarm in a turbulent flow in a
vertical square duct was addressed, and comparisons between unladen and bubble
laden flows were conducted. A DNS of unladen flow for a frictional Reynolds
number of $360$ was conducted, and unladen flow turbulence quantities were
validated with the literature. The influence of bubbles on liquid phase was
investigated, and the following overall conclusions can be drawn.

  * •

The introduction of bubbles strongly modifies the elongated flow structure
usually observed in the unladen flow, and the longest structures are seen to
occupy half of the duct. In the laden flow, instantaneous flow structures are
oscillatory, and flatter. The longer structures are concentrated near the
wall, and the shorter ones are distributed randomly away from the wall. The
instantaneous laden flow is dominated by high and low values of positive and
negative fluctuations respectively.

  * •

Starting from a uniform configuration, bubbles migrate toward the wall. They
form a layer near the wall and a cluster near the corner. The lateral
migration of bubbles, which leads to preferential concentration, is a result
of balance between shear-induced lift force and turbophoresis effects. The
secondary flows, which transport high energy fluids from the central region to
the wall region along the corner bisector, are also responsible for the
migration of bubbles in the corner region.

  * •

Turbulence quantities, including velocity fluctuations and Reynolds stress,
are enhanced by the introduction of bubbles. In the unladen flow, all
turbulence quantities have one peak and one valley. Along the wall bisector
they achieved their peaks near the wall and valleys at the center, but in the
laden flow, multiple peaks and valleys were observed. The first peak was
observed near the wall, and the second peak was observed between the wall and
half distance.

## Acknowledgment

The authors thank the financial support from the Air Conditioning and
Refrigeration Center (ACRC), University of Illinois at Urbana-Champaign, USA.
This research is also part of the Blue Waters sustained-petascale computing
project, which is supported by the National Science Foundation (awards
OCI-0725070 and ACI-1238993) and the state of Illinois. Blue Waters is a joint
effort of the University of Illinois at Urbana-Champaign and its National
Center for Supercomputing Applications. The authors further thank NVIDIA
Hardware Grant Program for providing the GPUs for an in-house workstation.

## References

  * Tomiyama et al. [2002] A. Tomiyama, H. Tamai, I. Zun, and S. Hosokawa.  Transverse migration of single bubbles in simple shear flows.  _Chemical Engineering Science_ , 57:1849–1858, 2002.  ISSN 0066-4189.  doi: 10.1146/annurev.fluid.37.061903.175803.
  * Afghan [1988] N. H. Afghan, editor.  _Transient Phenomena in Multiphase Flow_.  Hemisphere Publishing, 1988.
  * Zun et al. [1993] I. Zun, I. Kljenak, A. Serizawa, and S. Moze.  Space time evolution of non homogeneous bubble distribution in upward flow.  _International Journal of Multiphase Flow_ , 19:151–172, 1993.  ISSN 0066-4189.  doi: 10.1146/annurev.fluid.37.061903.175803.
  * Liu [1989] T.J. Liu.  _Experimental investigation of turbulence structure in two-phase bubbly flow_.  Phd thesis, Northwest University, Evanston, Illinois, 1989.
  * Liu and Bankoff [1993a] T. J. Liu and S. G. Bankoff.  Structure of air-water bubbly flow in a vertical pipe-1. Liquid velocity and turbulence measurements.  _International Journal of Heat and Mass Transfer_ , 36:1049–1060, 1993a.  ISSN 0066-4189.  doi: 10.1146/annurev.fluid.37.061903.175803.
  * Liu and Bankoff [1993b] T. J. Liu and S. G. Bankoff.  Structure of air-water bubbly flow in a vertical pipe-1. Void fraction, bubble velocity and bubble size distribution.  _International Journal of Heat and Mass Transfer_ , 36:1061–1072, 1993b.  ISSN 0066-4189.  doi: 10.1146/annurev.fluid.37.061903.175803.
  * Colin et al. [2012] Catherine Colin, Jean Fabre, and Arjan Kamp.  Turbulent bubbly flow in pipe under gravity and microgravity conditions.  _Journal of Fluid Mechanics_ , 711:469–515, 2012.
  * Lance and Bataille [1991] M. Lance and J. Bataille.  Turbulence in the liquid phase of a uniform bubbly air–water flow.  _Journal of Fluid Mechanics_ , 222:95–118, 1991.  ISSN 0022-1120.  doi: 10.1017/S0022112091001015.
  * Yuan and Michaelides [1992] Z. Yuan and E. E. Michaelides.  Turbulence modulation in particulate flows-A theoretical approach.  _International Journal of Multiphase Flow_ , 18(5):779–785, 1992.  ISSN 03019322.  doi: 10.1016/0301-9322(92)90045-I.
  * Hosokawa and Tomiyama [2013] Shigeo Hosokawa and Akio Tomiyama.  Bubble-induced pseudo turbulence in laminar pipe flows.  _International Journal of Heat and Fluid Flow_ , 40:97–105, 2013.
  * Alméras et al. [2015] Elise Alméras, Frédéric Risso, Véronique Roig, Sébastien Cazin, Cécile Plais, and Frédéric Augier.  Mixing by bubble-induced turbulence.  _Journal of Fluid Mechanics_ , 776:458–474, 2015.
  * Hirt and Nichols [1981] C. W Hirt and B. D Nichols.  Volume of Fluid (VOF) method for the dynamics of free boundaries.  _Journal of Computational Physics_ , 39(1):201–225, January 1981.
  * Rudman [1998] Murray Rudman.  A volume-tracking method for incompressible multifluid flows with large density variations.  _International Journal for Numerical Methods in Fluids_ , 28(2):357–378, August 1998.  ISSN 0271-2091.
  * Cummins et al. [2005] Sharen J. Cummins, Marianne M. Francois, and Douglas B. Kothe.  Estimating curvature from volume fractions.  _Computers & Structures_, 83(6-7):425 – 434, 2005.  ISSN 0045-7949.  doi: http://dx.doi.org/10.1016/j.compstruc.2004.08.017.  URL http://www.sciencedirect.com/science/article/pii/S0045794904004110.
  * Francois et al. [2006] Marianne M. Francois, Sharen J. Cummins, Edward D. Dendy, Douglas B. Kothe, James M. Sicilian, and Matthew W. Williams.  A balanced-force algorithm for continuous and sharp interfacial surface tension models within a volume tracking framework.  _Journal of Computational Physics_ , 213(1):141–173, March 2006.  ISSN 00219991.
  * Wang and Tong [2010] Zhaoyuan Wang and Albert Y. Tong.  A sharp surface tension modeling method for two-phase incompressible interfacial flows.  _International Journal for Numerical Methods in Fluids_ , 64(7):709–732, September 2010.  ISSN 02712091.
  * Kumar and Vanka [2015] P. Kumar and S.P. Vanka.  Effects of confinement on bubble dynamics in a square duct.  _International Journal of Multiphase Flow_ , 77:32 – 47, 2015.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.06.014.
  * Kumar et al. [2015a] Purushotam Kumar, Kai Jin, and S. Pratap Vanka.  A multi-GPU based accurate algorithm for simulations of gas-liquid flows.  _Computational Methods and Tools in Thermal Fluids Sciences_ , pages 1–17, 2015a.  doi: 10.1615/TFESC1.cmd.012709.  URL http://dl.astfe.org/conferences/tfesc,5451756e07e31e0f,768908a902b726f9.html.
  * Jin et al. [2016] K. Jin, P. Kumar, S. P. Vanka, and B. G. Thomas.  Rise of an argon bubble in liquid steel in the presence of a transverse magnetic field.  _Physics of Fluids_ , 28(9), 2016.  ISSN 10897666.  doi: 10.1063/1.4961561.
  * Kumar et al. [2015b] Purushotam Kumar, Kai Jin, and S. Pratap Vanka.  Bubble rise and deformation in a non-newtonian fluid in a square duct.  _Computational Methods and Tools in Thermal Fluids Sciences_ , pages 1–15, 2015b.  doi: 10.1615/TFESC1.cmd.012708.  URL http://dl.astfe.org/conferences/tfesc,5451756e07e31e0f,768908a902b726f9.html.
  * Kumar et al. [2019] Purushotam Kumar, Kai Jin, and Surya Pratap Vanka.  Numerical Simulation of a Gas Bubble Rising in Power-Law Fluids Using a Sharp Surface Force Implementation.  Proceedings of the ASME-JSME-KSME 2019 8th Joint Fluids Engineering Conference. Volume 5: Multiphase Flow. San Francisco, California, USA, July 28-August 01, 2019.  doi: 10.1115/AJKFluids2019-4769.
  * Rider and Kothe [1998] William J. Rider and Douglas B. Kothe.  Reconstructing volume tracking.  _Journal of Computational Physics_ , 141(2):112–152, April 1998.  ISSN 00219991.
  * Noh and Woodward [1976] WF Noh and P Woodward.  SLIC (simple line interface calculation).  In _Proceedings of the Fifth International Conference_ , pages 330–340, Berlin, 1976. Springer-Verlag.
  * Li [1995] Jie Li.  Calcul d’Interface Affine par Morceaux(piecewise linear interface calculation).  _C. R. Acad. Sci. Paris_ , 320(8):391–396, 1995\.
  * Ashgriz and Poo [1991] N Ashgriz and JY Poo.  FLAIR: flux line-segment model for advection and interface reconstruction.  _Journal of Computational Physics_ , 93(2):449–468, April 1991.  ISSN 00219991.
  * Renardy and Renardy [2002] Yuriko Renardy and Michael Renardy.  PROST: a parabolic reconstruction of surface tension for the volume-of-fluid method.  _Journal of Computational Physics_ , 183(2):400–421, December 2002.  ISSN 00219991.
  * Sussman [2003] Mark Sussman.  A second order coupled level set and volume-of-fluid method for computing growth and collapse of vapor bubbles.  _Journal of Computational Physics_ , 187(1):110–136, May 2003.  ISSN 00219991.
  * Sussman et al. [2007] M. Sussman, K.M. Smith, M.Y. Hussaini, M. Ohta, and R. Zhi-Wei.  A sharp interface method for incompressible two-phase flows.  _Journal of Computational Physics_ , 221(2):469–505, February 2007.  ISSN 00219991.
  * Gueyffier et al. [1999] Denis Gueyffier, Jie Li, Ali Nadim, Ruben Scardovelli, and Stéphane Zaleski.  Volume-of-Fluid interface tracking with smoothed surface stress methods for three-dimensional flows.  _Journal of Computational Physics_ , 152(2):423–456, July 1999.  ISSN 00219991.
  * Nickolls and Dally [2010] J. Nickolls and W. J. Dally.  The gpu computing era.  _IEEE Micro_ , 30(2):56–69, March 2010.  ISSN 0272-1732.  doi: 10.1109/MM.2010.41.
  * Cud [2015] CUDA programming guide, September 2015.  URL http://docs.nvidia.com/cuda/pdf/CUDA_C_Programming_Guide.pdf.
  * Nickolls et al. [2008] John Nickolls, Ian Buck, Michael Garland, and Kevin Skadron.  Scalable parallel programming with cuda.  _Queue_ , 6(2):40–53, March 2008.  ISSN 1542-7730.  doi: 10.1145/1365490.1365500.  URL http://doi.acm.org/10.1145/1365490.1365500.
  * Association [2000] I. T. Association.  _InfiniBand Architecture Specification: Release 1.0_.  InfiniBand Trade Association, 2000.
  * Gropp et al. [1996] W Gropp, E Lusk, N Doss, and A Skjellum.  A high-performance, portable implementation of the mpi message passing interface standard.  _Parallel Computing_ , 22(6):789–828, 1996.
  * Gropp et al. [1999] W Gropp, E Lusk, and R Thakur.  _Using MPI-2: Advanced Features of the Message Passing Interface_.  MIT press, 1999.
  * Kumar [2016] Purushotam Kumar.  _Development and application of Lattice Boltzmann and accurate volume of fluid numerical techniques on graphics processing units_.  PhD thesis, University of Illinois at Urbana-Champaign, 2016.
  * Horwitz et al. [2014] J.A.K. Horwitz, P. Kumar, and S.P. Vanka.  Three-dimensional deformation of a spherical droplet in a square duct flow at moderate reynolds numbers.  _International Journal of Multiphase Flow_ , 67:10 – 24, 2014.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/j.ijmultiphaseflow.2014.07.010.  URL http://www.sciencedirect.com/science/article/pii/S0301932214001396.
  * Horwitz et al. [2019] Jeremy A. K. Horwitz, S. P. Vanka, and P. Kumar.  Lbm simulations of dispersed multiphase flows in a channel: Role of a pressure poisson equation.  Proceedings of the ASME-JSME-KSME 2019 8th Joint Fluids Engineering Conference. Volume 5: Multiphase Flow. San Francisco, California, USA, July 28-August 01, 2019.  doi: 10.1115/AJKFluids2019-4943.
  * Kim et al. [1987] John Kim, Parviz Moin, and Robert Moser.  Turbulence statistics in fully developed channel flow at low reynolds number.  _Journal of Fluid Mechanics_ , 177:133–166, 4 1987.  ISSN 1469-7645.  doi: 10.1017/S0022112087000892.  URL http://journals.cambridge.org/article_S0022112087000892.
  * Madabhushi and Vanka [1991] Ravi K Madabhushi and S P Vanka.  Large eddy simulation of turbulence-driven secondary flow in a square duct.  _Physics of Fluids A: Fluid Dynamics_ , 3(11):2734, 1991.  ISSN 08998213.  doi: 10.1063/1.858163.
  * Panidis and Papailiou [2000] Thrasyvoulos Panidis and Demosthenes D. Papailiou.  The structure of two-phase grid turbulence in a rectangular channel: An experimental study.  _International Journal of Multiphase Flow_ , 26(8):1369–1400, 2000.  ISSN 03019322.  doi: 10.1016/S0301-9322(99)00085-3.
  * Santarelli and Fröhlich [2015] C. Santarelli and J. Fröhlich.  Direct Numerical Simulations of spherical bubbles in vertical turbulent channel flow.  _International Journal of Multiphase Flow_ , 75:174–193, 2015.
  * Bunner and Tryggvason [2002] B Bunner and G. Tryggvason.  Dynamics of homogeneous bubbly flows Part 2. Velocity fluctuations.  _Journal of Fluid Mechanics_ , 466:53–84, 2002.  ISSN 0022-1120.  doi: 10.1017/S0022112002001180.
  * Bunner and Tryggvason [2003] B. Bunner and G. Tryggvason.  Effect of bubble deformation on the properties of bubbly flows.  _Journal of Fluid Mechanics_ , 495:77–118, 2003.  ISSN 00221120.  doi: 10.1017/S0022112003006293.
  * Grace [1973] JR Grace.  Shapes and velocities of bubbles rising in infinite liquids.  _Trans. Inst. Chem. Eng._ , 51(2):116–120, 1973\.
  * Gavrilakis [1992] S. Gavrilakis.  Numerical simulation of low-Reynolds-number turbulent flow through a straight square duct.  _Journal of Fluid Mechanics_ , 244(-1):101, 1992\.  ISSN 0022-1120.  doi: 10.1017/S0022112092002982.
  * Huser and Biringen [1993] Asmund Huser and Sedat Biringen.  Direct numerical simulation of turbulent flow in a square duct.  _Journal of Fluid Mechanics_ , 257:65–95, 12 1993.  ISSN 1469-7645.  doi: 10.1017/S002211209300299X.  URL http://journals.cambridge.org/article_S002211209300299X.
  * Broglia et al. [2003] Riccardo Broglia, Andrea Pascarelli, and Ugo Piomelli.  Large-eddy simulations of ducts with a free surface.  _Journal of Fluid Mechanics_ , 484:223–253, 2003.  ISSN 00221120.  doi: 10.1017/S0022112003004257.
  * Zhang et al. [2015] Hao Zhang, F. Xavier Trias, Andrey Gorobets, Yuanqiang Tan, and Assensi Oliva.  Direct numerical simulation of a fully developed turbulent square duct flow up to $Re_{\tau}$ = 1200.  _International Journal of Heat and Fluid Flow_ , 54:258–267, 2015.
  * Winkler et al. [2004] C.M. Winkler, Sarma L. Rani, and S.P. Vanka.  Preferential concentration of particles in a fully developed turbulent square duct flow.  _International Journal of Multiphase Flow_ , 30(1):27 – 50, 2004.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/j.ijmultiphaseflow.2003.11.003.  URL http://www.sciencedirect.com/science/article/pii/S0301932203002039.
  * Winkler et al. [2006] C.M. Winkler, Sarma L. Rani, and S.P. Vanka.  A numerical study of particle wall-deposition in a turbulent square duct flow.  _Powder Technology_ , 170(1):12 – 25, 2006.  ISSN 0032-5910.  doi: http://dx.doi.org/10.1016/j.powtec.2006.08.009.  URL http://www.sciencedirect.com/science/article/pii/S0032591006003391.
  * Chaudhary et al. [2009] R. Chaudhary, S. P. Vanka, and B. G. Thomas.  Direct numerical simulations of magnetic field effects on turbulent duct flows.  In _Proceedings of the ASME 2009 International Mechanical Engineering Congress and Exposition (IMECE 2009)_ , pages 229–308, 2009.  doi: 10.1115/IMECE2009-11527.  URL http://proceedings.asmedigitalcollection.asme.org/proceeding.aspx?articleid=1640392.
  * Chaudhary et al. [2010] R. Chaudhary, S. P. Vanka, and B. G. Thomas.  Direct numerical simulations of magnetic field effects on turbulent flow in a square duct.  _Physics of Fluids_ , 22(7):075102, 2010.  doi: http://dx.doi.org/10.1063/1.3456724.  URL http://scitation.aip.org/content/aip/journal/pof2/22/7/10.1063/1.3456724.
  * Chaudhary et al. [2011] R. Chaudhary, A.F. Shinn, S.P. Vanka, and B.G. Thomas.  Direct numerical simulations of transverse and spanwise magnetic field effects on turbulent flow in a 2:1 aspect ratio rectangular duct.  _Computers & Fluids_, 51(1):100 – 114, 2011\.  ISSN 0045-7930.  doi: http://dx.doi.org/10.1016/j.compfluid.2011.08.002.  URL http://www.sciencedirect.com/science/article/pii/S004579301100243X.
  * Shinn [2011] Aaron F. Shinn.  _Large eddy simulations of turbulent flows on graphics processing units: application to film-cooling flows_.  PhD thesis, University of Illinois at Urbana-Champaign, 2011.