.tifpng.pngconvert #1 .tif

# Dynamics of two air bubbles rising in a shear-thinning fluid

Purushotam Kumar<EMAIL_ADDRESS>Corning Incorporated, Manufacturing
Technology & Engineering, Corning, NY 14830  Surya P. Vanka
<EMAIL_ADDRESS>Mechanical Science & Engineering, University of Illinois
Urbana-Champaign, Urbana, IL 61801

###### Abstract

In this paper, we have studied the three-dimensional dynamics of two equally
sized air bubbles rising in a shear-thinning fluid. We have used the combined
level set and volume of fluid (CLSVOF) method to track interface, maintain
mass balance and estimate the interface curvature. Additionally, we have
incorporated a Sharp Surface Force Method (SSF) for surface tension forces.
This method significantly suppressed the spurious velocities commonly observed
with the conventional volume of fluid method and the Continuum Surface Force
(CSF). The algorithm is implemented in an in-house code called CUFLOW and runs
on multiple GPUs platform.

We have explored the effects of fluid rheology on the three-dimensional
dynamics of two in-line bubbles. Two power-law indices (0.5 and 1) are
investigated to highlight differences in shear-thinning and Newtonian fluids.
For a range of parameters examined here, bubbles motion in a shear-thinning
fluid is seen to be unsteady with significant shape oscillations. Further, we
have examined the rise velocity, droplets rise path, transient shapes and
found that modification of viscosity by the motion of the leading bubble
changes the dynamics of the trailing bubble.

Keywords: Bubble dynamics; Bubble interactions; Liquid flow structure; Two-
phase flows; Volume of fluid (VOF) method

##  1 Introduction

A gas bubble rising in a liquid medium is one of the simplest two-phase flow
problems and addresses several important interactions between two phases. This
problem has been previously studied to understand lift Morsi and Alexander
(1972); Kurose et al. (2001) and drag Magnaudet and Legendre (1998) forces
acting on a bubble, lateral migration of the bubble in a co-flow Adoua et al.
(2009), and bubble-wall and bubble-wake interactions. A large number of
researchers, including Allen (1900); Hadamard (1911); Rybczynski (1911); Bond
(1927); Moore (1959) have investigated the bubble deformation, rise path,
terminal velocity, and confinement effects through experiments, theoretical
modeling, and numerical simulations.

Dynamics of a rising bubble are a function of the Eötvös number
($Eo=\Delta\rho gd^{2}\sigma^{-1}$, also referred as Bond number
$Bo=\rho_{l}gd^{2}\sigma^{-1}$), Morton number
$(Mo=g\mu_{o}^{4}\rho_{l}^{-1}\sigma^{-3})$, confinement ratio
$(C_{r}=Wd^{-1})$, and other external forces such as a magnetic field Jin2016
and variable viscosity Kumar et al. (2015a). An exhaustive collection of
experimental data for a bubble rising in an unconfined medium has been
summarized in works of Grace (1973); Grace et al. (1976); Bhaga and Weber
(1981), available as nomogram charts of the terminal Reynolds and bubble shape
for combinations of Bond and Morton numbers. In addition, the proximity of the
bubble to a wall alters the bubble rise velocity and shape and increases the
shear force at the gas-liquid interface Figueroa-Espinoza et al. (2008); Kumar
and Vanka (2015). Detailed reviews of bubble dynamics in Newtonian fluids are
available in the works of Clift et al. (1978); DeKee (2002).

The dynamics of bubbles in Newtonian fluids have been widely studied. However,
many fluids have non-Newtonian behavior and are used in the chemical and
pharmaceutical industries. For instance, the fermentation of liquor, sugarcane
processing, ethanol extraction, aeration of the membrane is typical shear-
thinning applications. Bubble rise in non-Newtonian fluids is more complex
because the bubble-generated velocities subsequently change the local fluid
viscosity due to the shear. Further, the nature of the fluid (Bingham,
thixotropic, shear-thinning, shear-thickening, etc.) is important in modifying
the fluid viscosity in the bubble vicinity and thereby the local shear stress
values. Thus, a large parameter space is available, even for a single bubble
in a non-Newtonian fluid, which provides rich and complex flow physics
regions. Previously, several efforts have been made towards understanding such
bubble motion in viscoelastic and inelastic non-Newtonian fluids. Some example
works are by Astarita and Apuzzo (1965); Leal et al. (1971); Carreau et al.
(1974); Zana and Leal (1978); Kawase and Ulbrecht (1981); Dekée et al. (1986);
Gummalam et al. (1988); Kee et al. (1990); Manjunath and Chhabra (1992);
Miyahara and Yamanaka (1993); Liu et al. (1995); Rodrigue et al. (1996); Li et
al. (1997); Chhabra (1998); Gauri and Koelling (1998); Dubash and Frigaard
(2007); Tsamopoulos et al. (2008); Sikorski et al. (2009); Kishore et al.
(2013); Tripathi et al. (2015) for bubble motion in shear-thinning, shear-
thickening and Bingham fluids. The papers of Kulkarni and Joshi (2005);
Chhabra (2006) provide good summaries of the important results.

The dynamics of multiple bubbles rising has been studied by several authors
(Abbassi et al., 2018; Gumulya et al., 2017; Watanabe and Sanada, 2006; Yuan
and Prosperetti, 1994; Legendre et al., 2003; Yu et al., 2011; Katz and
Meneveau, 1996; Ruzicka, 2000) and have reported that flow structure around
the bubble has an important impact on the interaction between leading and
trailing bubbles. They have reported that wake in front of a trailing bubble
plays an important role in reducing the drag around it. Due to reduced drag,
the trailing bubble accelerates until it coalesces with the leading bubble.
Along with these, there are other studies; however, the majority of the
studies for multiple bubble rise are conducted in the Newtonian fluids.

In the present paper, we have briefly discussed a numerical technique to
simulate two-phase flows and used it to study the three-dimensional
deformation and rise of two equally sized bubbles in shear-thinning and
Newtonian fluids for a given bubble size and center to center distance.
Despite a reasonably large number of previous studies, there are very few
studies that have examined the dynamics of in-line bubbles rising in non-
Newtonian fluids. Various quantities such as the bubble shape, rise velocity,
and rise paths are investigated. The algorithm used for the computations is
presented in the numerical method section. The problem description is provided
in the computational details section 3. In the results and discussion section
4, we present the results of the simulations and discuss the critical
findings. A summary of the present findings is given at the end.

##  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, the inclusion of surface tension forces,
and the solution of Navier-Stokes equations. The volume of fluid (VOF) method
of Hirt and Nichols Hirt and Nichols (1981) is used for interface capturing.
The truncated interface between the gas and liquid in a control volume is
represented as an oblique plane, and the interface normal, and curvature are
computed using a smoothed liquid volume fraction and the height function
method (Rudman Rudman (1998) and Cummins et al. Cummins et al. (2005)),
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.Francois et al. (2006), Wang and Tong Wang and Tong (2010)).
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 Kumar and Vanka
(2015); Kumar et al. (2015b) and the effects of magnetic field on bubble
dynamics Jin et al. (2016).

### Governing Equations

We assume that both the gas and liquid are isothermal and incompressible. 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
the 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 the fluid velocity, $p$ is pressure,
$\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 the acceleration due to gravity.

Since the 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 Rider and Kothe
(1998)) coupled with a second-order operator split method (Noh and Woodward
Noh and Woodward (1976), Li Li (1995), Ashgriz and Poo Ashgriz and Poo (1991)
and Francois et al. Francois et al. (2006)) is used to solve the evolution
equation. The mixture density and viscosity are calculated as a 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 reduces 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 Renardy and
Renardy (2002), Sussman et al. Sussman (2003); Sussman et al. (2007),
Gueyffier et al. Gueyffier et al. (1999)); however, in our experience, the
current method is more accurate for reduction of spurious velocities and
relatively easier to implement on GPUs.

### 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}=$
$\displaystyle\rho_{c}^{n}\mathbf{u}_{c}^{n}-\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 are described in Kumar and
Vanka Kumar and Vanka (2015).

The entire algorithm has been programmed on the CUDA-Fortran platform
supported by the Portland Group (PGI) Fortran compiler to run on multiple
graphics processing units (GPU). Previously, we have validated the multiple
GPU implementation (Kumar et al., 2015b; Vanka et al., 2016; Kumar, 2016) to
study confinement effects on bubble dynamics (Kumar and Vanka, 2015), two-
phase flows at T-juctions (Horwitz et al., 2012, 2013; Kumar et al., 2013),
bubble dynamics in non-Newtonain fluids (Kumar et al., 2015a, 2019), droplet
dynamics in square duct (Horwitz et al., 2014, 2019), Argon bubble rising in
liquid steel (Vanka et al., 2015; Jin et al., 2016; Kumar and Vanka, 2022,
2023) and turbulent bubbly flow (Vanka et al., 2016; Kumar and Vanka, 2021).

##  3 Computational Details

### Problem Setup

Figure 1 shows the geometry of the problem considered. The dimensions of
domain are $4d\times 4d\times 24d$ along $x$, $y$ and $z$ directions
respectively. The height was chosen such that the bubble in a Newtonian fluid
can reach a steady velocity at the end of the domain. All boundaries of the
domain are considered to be no-slip and no-penetration walls. The trailing
bubble is initially placed two diameters above the bottom boundary in the
center of the duct cross-section. The leading bubble is placed according to
the bubble’s center to center distance of the problem (in this case $4d$). The
gravitational force acts downwards along the negative z-axis. Based on a
systematic grid refinement study, we have selected $32$ control volumes per
bubble diameter as an adequate resolution, with a grid of $128\times 128\times
768$ ($\approx$ 13 million) control volumes current study.

Figure 1: Initial position of the bubble

The liquid and gas densities are taken as 1000 and 1.2 $kg/m^{3}$. The zero
shear viscosity $(\mu_{0})$ and time constant $(\lambda)$ of the liquid are
$0.014$ $Pa.s$ and 1 $s$ respectively. The power-law index is varied from each
case. The dynamic viscosity of the gas $(\mu_{g})$ is taken as $1.8\times
10^{-5}Pa.s$. The surface tension $(\sigma)$ between liquid and gas phase is
prescribed a value of $0.072N/m$111In practice there is variation in the
surface tension due to addition of polymers for preparation of shear-thinning
and shear-thickening fluids but the variation is relatively less significant
for low concentration of polymers Ohta et al. (2003); Zhang et al. (2010);
Ohta et al. (2012).

The liquid viscosity is calculated using the Carreau model Carreau (1972) with
the power-law index 0.5 and 1.0 for shear-thinning and constant viscosity,
respectively. It is calculated using,

$\displaystyle\mu_{l}(\dot{\gamma})=\mu_{o}(1+\left(\lambda\dot{\gamma}\right)^{2})^{\nicefrac{{n-1}}{{2}}}$
(8)

Where, $\mu_{o}$, $\lambda$ and $n$ are zero shear viscosity, time constant
and power-law index, respectively. The shear rate $(\dot{\gamma})$ is given
as,

$\displaystyle\dot{\gamma}=\sqrt{2\mathbf{D}:\mathbf{D}}$ (9)

The deformation tensor $\mathbf{D}$ is defined as,

$\displaystyle\mathbf{D}=\nicefrac{{1}}{{2}}\left(\nabla\textbf{u}+\nabla\textbf{u}^{T}\right)$
(10)

The dynamics of two bubbles ascending are governed by four non-dimensional
parameters,

Eötvös number $\displaystyle=Eo=\frac{\Delta\rho gd^{2}}{\sigma}$ (11) Morton
number
$\displaystyle=Mo=\frac{g\mu_{l}^{4}\Delta\rho}{\rho_{l}^{2}\sigma^{3}}$ (12)
Confinement ratio $\displaystyle=C_{r}=\frac{W}{d}$ (13) center to center
distance $\displaystyle=h_{d}=\frac{h}{d}$ (14)

In the above equations, $\Delta\rho=\rho_{l}-\rho_{g}$, $W$ is width of liquid
column and $h$ is the initial center to center distance of bubbles. For a gas-
liquid system the density ratio is large, hence
$\Delta\rho=\rho_{l}-\rho_{g}\approx\rho_{l}$ and Eötvös number can be also
represented by Bond number. In the current study, the viscosity of the liquid
phase $(\mu_{l})$ varies with the shear rate $(\dot{\gamma})$, therefore we
define the Morton number using the zero-shear liquid viscosity. Hence, the key
non-dimensional parameters are given by:

Bond number $\displaystyle=Bo=\frac{\rho_{l}gd^{2}}{\sigma}$ (15) Morton
number $\displaystyle=Mo=\frac{g\mu_{l,0}^{4}}{\rho_{l}\sigma^{3}}$ (16)
Confinement ratio $\displaystyle=C_{r}=\frac{W}{d}$ (17) center to center
distance $\displaystyle=h_{d}=\frac{h}{d}$ (18)

##  4 Results and Discussion

We now present the results of a study to analyze the effects of fluid rheology
on bubble deformation, rise velocity, rise path, and liquid viscosity
distribution. Two power-law indices (0.5 and 1.0) have been considered to
understand its effects. The power-law indices capture the shear-thinning and
Newtonian behavior of the non-Newtonian fluid. The zero-shear viscosity
$(\mu_{0})$ and time constant $(\lambda)$ of the Carreau model Carreau (1972)
have been kept constant to isolate the effects of the power-law index. The
bubbles are initially placed at the center of the cross-section, but the
bubbles may come closer to the wall during the ascend. The Bond is kept
constant at $Bo=2$, and the terminal shape in the Newtonian fluid is expected
to be ellipsoidal.

Now, we discuss the effects of the power-law index $(n)$ for a center to
center distance of $4d$. For two liquids considered here, the zero-shear
Morton number $(Mo)$ is $10^{-6}$. According to the Grace diagram Grace
(1973), an initially spherical bubble of $Bo=2$ rising in a Newtonian fluid
with a Morton number of $10^{-6}$ deforms into an ellipsoidal shape. In non-
Newtonian fluids, the viscosity will be modified by the shear generated
adjacent to the bubble. Hence, the effective viscosity will be lower for the
shear-thinning fluid. Consequently, the effective Morton number of shear-
thinning fluid will be lower than that of the Newtonian fluid.

|  | time (ms)
---|---|---
$n$ |  |  0 |  80 |  160 |  240 |  320
0.5 |  LB |  |  |  |  |
0.5 |  TB |  |  |  |  |
1.0 |  LB |  |  |  |  |
1.0 |  TB |  |  |  |  |
Figure 2: Transient three-dimensional front views of the bubbles for center to
center distance = $4d$, and two power-law indices $(n)$

Figure 2 shows the transient shapes of leading (LB) and trailing (TB) bubbles
for $n$ = 0.5 and 1 and $h_{d}$ = 4. From fig. 2 we can observe that the
bubble terminal shape in the Newtonian fluid (Fig. 2 row 3) is ellipsoidal
with nearly symmetrical top and bottom surfaces. However, in the shear-
thinning fluid (Fig. 2 row 1), the bubble deforms into an ellipsoid that is
flatter than its Newtonian counterpart with definite asymmetry between top and
bottom surfaces. It can also be noticed that the bubble in shear-thinning
fluid goes through shape and orientation oscillations. In contrast, in the
Newtonian fluid, the bubble achieves a steady shape and rises without any
orientation changes.

For the simulation time considered in this study, the dynamics of leading and
trailing bubbles are quite similar in a Newtonian fluid. On the other hand, in
shear-thinning fluid, their dynamics are complex. At $t=80$ $ms$, the shape of
both bubbles is nearly identical, primarily because the trailing bubble has
not yet reached the flow field modified by the leading bubble. The orientation
and shape of the leading and trailing bubbles are different afterward. At
$t=160$ $ms$, the surface normal at the top is upward-front and upward-back
orientated for leading and trailing bubbles, respectively. This indicates that
vortex shedding behind the bubble is oriented differently. At $t=320$ $ms$,
the shape of the trailing bubble is upside-down of the leading bubble.

(a) $n$ = 0.5

(b) $n$ = 1.0

Figure 3: Bubble rise path for center to center distance = $4d$, and two
power-law indices

Figure 3 shows the ascending path of leading and trailing bubbles for a center
to center distance of $4d$ and two power-law indices. It can be noticed that
bubbles in the Newtonian fluid rise in a rectilinear path. This is in line
with our observation from Fig. 2, the bubble reached a steady shape without
any shape or orientation oscillations. It has also been reported earlier by
Kumar et al. [???] that in the absence of bubble shape or orientation
oscillations, the vortices behind the bubble are symmetric and do not amplify
any instability. We are also observing similar phenomena for both bubbles in
the Newtonian fluid.

In the shear-thinning fluid, the rise path of bubbles goes from rectilinear to
oscillatory. It should be noted that initially the trailing bubble has a
smaller velocity and is away from the flow field modified by the leading
bubble. Its shape changes from the spherical to asymmetric ellipsoid and
doesn’t exhibit surface oscillations. This is the primary reason for the
initial rectilinear path of the bubble. As the bubble accelerates, it reaches
a certain critical speed (or Reynold number), after which it sheds asymmetric
vortices. This leads to non-uniform pressure distribution behind the bubble,
which pushes the bubble away from the center of the liquid column and creates
a zig-zag rise path. From Fig. 3(a) we can observe that the leading bubble
initially migrates towards the left compared to the trailing bubble migrating
towards the right.

Figure 4: Rise velocity of bubbles for center to center distance = $4d$, and
two power-law indices

Figure 4 shows the rise velocities of leading and trailing bubbles in the
Newtonian and shear-thinning fluids for a center to center distance of 4. In
the Newtonian fluid, the velocity curves of both bubbles are identical for up
to $t\approx 60\leavevmode\nobreak\ ms$, and afterward, the curves diverge.
The velocity of the leading bubble between $60\leavevmode\nobreak\ ms$ and
$200\leavevmode\nobreak\ ms$ is almost constant, whereas the trailing bubble
rapidly accelerates to achieve a higher velocity than the leading bubble. This
behavior is primarily because the leading bubble ascends through the
undisturbed flow field hence experiences maximum drag. The trailing bubble
ascends through a modified flow field; as a result, it experiences lower drag
force and achieves higher rise velocity.

Similar to the Newtonian case, rise velocity curves in the shear-thinning
fluid are identical for up to $t\approx 70\leavevmode\nobreak\ ms$, and
afterward, the curves diverge. The rise velocity goes through rapid
oscillations. In the Newtonian fluid, the trailing bubble velocity was always
higher than that of the leading bubble. However, in the shear-thinning fluid,
there are instances where the leading bubble rises with higher velocity than
the trailing bubble. It can be noted for any given time window, the amplitude
of oscillation of trailing bubble rise velocity is higher than that of the
leading bubble.

(a) $n=0.5$

(b) $n=1.0$

Figure 5: Velocity vector at $t=160\leavevmode\nobreak\ ms$ for center to
center distance = $4d$, and two power-law indices

Figures 5(a) and 5(b) show the velocity vectors on $yz-$plane at
$t=160\leavevmode\nobreak\ ms$ for shear-thinning and constant viscosity
cases. The interface has been drawn to show the cross-section of the bubble on
the $yz-$ plane. We notice from Fig. 5(b) that velocity vectors in front of
the trailing bubble are streamlined without any recirculation zones, and on
both left and right sides, the flow distribution is symmetric. From Fig. 5(a)
we can notice that the flow field in front of the trailing bubble contains
recirculation zones and is asymmetric on the left and right sides of the
bubble. The interaction between the recirculation zone and bubble creates
asymmetric forces and pushes away from the centerline. The impact of this was
observed in Fig. 3(a).

(a) Leading bubble

(b) Trailing bubble

Figure 6: Velocity vector at $t=160\leavevmode\nobreak\ ms$ for center to
center distance = $4d$, and two power-law index = 0.5

(a) Leading bubble

(b) Trailing bubble

Figure 7: Velocity vector at $t=160\leavevmode\nobreak\ ms$ for center to
center distance = $4d$, and two power-law index = 1.0

Figures 6 and 7 show the zoomed-in view of bubbles at
$t=160\leavevmode\nobreak\ ms$ on the $yz-$plane. The zoom ratio of the four
subfigures is kept the same to show that the deformation of the bubble in the
shear-thinning fluid is substantially higher than that in the Newtonian fluid.
We can further observe the asymmetric and symmetric flow field around the
trailing bubbles in shear-thinning and Newtonian fluid, respectively.

##  5 Conclusions

In this numerical study, we have compared the effects of fluid rheology on the
three-dimensional dynamics of two in-line bubbles rising in a quiescent liquid
medium. We have used a volume of fluid method coupled with a three-dimensional
geometry construction method for accurate representation and advection of the
interface. We have modeled the surface tension using the sharp surface force
(SSF) method. The fractional step method with the implicit treatment of the
viscous terms and geometry construction method for the convection term has
been used to solve the Navier-Stokes equations. The effects of the power-law
index on the bubble shapes, rise velocity, rise path, and flow field at
representative times have been explored.

We observe that the bubble deformation is strongly dependent on the power-law
index of the fluid – in the Newtonian fluid, both leading and trailing bubbles
deformed into an ellipsoid with symmetric top and bottom surfaces and attained
steady-state shape. In contrast, the bubble became more flatter with
asymmetric top and bottom surfaces in the shear-thinning fluid and went
through surface oscillations as they ascended.

Further, we observe that the bubbles initially rose in a linear path in shear-
thinning fluid and later followed a zigzag path. In contrast, in the Newtonian
fluid, they rose in a rectilinear path throughout the ascent. The zigzag
motion of the leading bubble was triggered by an asymmetry in the flow field
behind and consequent vortex shedding behind the bubble. The zigzag motion of
the trailing bubble is the combination of the unsteady flow in front of it and
the asymmetric flow around it.

The rise velocity oscillates in the shear-thinning fluids, whereas no such
features are present in the Newtonian and shear-thickening fluids. These
oscillations are related to the shape oscillation of the bubble and a zigzag
rise path. In both liquids, leading and trailing bubbles initially ascended
with identical speeds. Once they achieved maximum deformation, the trailing
bubble in the Newtonian fluid rose faster than the leading bubble. In the
shear-thinning fluid, the magnitude of oscillation of rising velocity was
higher for the trailing bubble than the leading bubble.

## References

  * Morsi and Alexander (1972) S. A. Morsi and A. J. Alexander.  An investigation of particle trajectories in two-phase flow systems.  _Journal of Fluid Mechanics_ , 55:193–208, 9 1972.  ISSN 1469-7645.  doi: 10.1017/S0022112072001806.
  * Kurose et al. (2001) R. Kurose, R. Misumi, and S. Komori.  Drag and lift forces acting on a spherical bubble in a linear shear flow.  _International Journal of Multiphase Flow_ , 27(7):1247 – 1258, 2001.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/S0301-9322(00)00073-2.
  * Magnaudet and Legendre (1998) Jacques Magnaudet and Dominique Legendre.  The viscous drag force on a spherical bubble with a time-dependent radius.  _Physics of Fluids_ , 10(3), 1998.
  * Adoua et al. (2009) Richard Adoua, Dominique Legendre, and Jacques Magnaudet.  Reversal of the lift force on an oblate bubble in a weakly viscous linear shear flow.  _Journal of Fluid Mechanics_ , 628:23–41, 6 2009.  ISSN 1469-7645.  doi: 10.1017/S0022112009006090.  URL http://journals.cambridge.org/article_S0022112009006090.
  * Allen (1900) H. S. Allen.  The motion of a sphere in a viscous fluid.  _Philosophical Magazine Series 5_ , 50(304):323–338, 1900.
  * Hadamard (1911) J. Hadamard.  Mouvement permanent lent d’une sphere liquide et visqueuse dans une liquide visqueux.  _Comptes rendus de l’Académie des sciences_ , 152:1735, 1911.
  * Rybczynski (1911) W. Rybczynski.  On the translatory motion of a fluid sphere in a viscous medium.  _Bull. Acad. Sci., Krakow, Series A_ , 40:40–46, 1911.
  * Bond (1927) W.N. Bond.  Bubbles and drops and Stokes’ law.  _Philosophical Magazine Series 7_ , 4(24):889–898, 1927.
  * Moore (1959) D. W. Moore.  The rise of a gas bubble in a viscous liquid.  _Journal of Fluid Mechanics_ , 6(01):113–130, March 1959.  ISSN 0022-1120.
  * Kumar et al. (2015a) 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, 2015a.
  * Grace (1973) JR Grace.  Shapes and velocities of bubbles rising in infinite liquids.  _Trans. Inst. Chem. Eng._ , 51(2):116–120, 1973\.
  * Grace et al. (1976) J.R. Grace, T. Wairegi, and T.H Nguyen.  Shapes and velocities of single drops and bubbles moving freely through immiscible liquids.  _Chemical Engineering Research and Design_ , 54a:167–173, 1976.
  * Bhaga and Weber (1981) Dahya Bhaga and M.E. Weber.  Bubbles in viscous liquids: shapes, wakes and velocities.  _Journal of Fluid Mechanics_ , 105:61–85, 1981.
  * Figueroa-Espinoza et al. (2008) B. Figueroa-Espinoza, R. Zenit, and D. Legendre.  The effect of confinement on the motion of a single clean bubble.  _Journal of Fluid Mechanics_ , 616:419–443, November 2008\.  ISSN 0022-1120.
  * 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.
  * Clift et al. (1978) Roland Clift, John R Grace, and Martin E Weber.  _Bubbles, Drops and Particles_.  Academic Press, London, 1978.
  * DeKee (2002) Daniel DeKee.  _Transport Processes in Bubbles, Drops and Particles_.  CRC Press, 2002.
  * Astarita and Apuzzo (1965) Gianni Astarita and Gennaro Apuzzo.  Motion of gas bubbles in non-newtonian liquids.  _AIChE Journal_ , 11(5):815–820, 1965.  ISSN 1547-5905.  doi: 10.1002/aic.690110514.  URL http://dx.doi.org/10.1002/aic.690110514.
  * Leal et al. (1971) L. G. Leal, J. Skoog, and A. Acrivos.  On the motion of gas bubbles in a viscoelastic liquid.  _The Canadian Journal of Chemical Engineering_ , 49(5):569–575, 1971.  ISSN 1939-019X.  doi: 10.1002/cjce.5450490504.  URL http://dx.doi.org/10.1002/cjce.5450490504.
  * Carreau et al. (1974) P. J. Carreau, M. Devic, and M. Kapellas.  Dynamique des bulles en milieu viscoélastique.  _Rheologica Acta_ , 1974.
  * Zana and Leal (1978) E. Zana and L.G. Leal.  The dynamics and dissolution of gas bubbles in a viscoelastic fluid.  _International Journal of Multiphase Flow_ , 4(3):237 – 262, 1978.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/0301-9322(78)90001-0.
  * Kawase and Ulbrecht (1981) Yoshinori Kawase and Jaromir J. Ulbrecht.  On the abrupt change of velocity of bubbles rising in non-newtonian liquids.  _Journal of Non-Newtonian Fluid Mechanics_ , 8(3–4):203 – 212, 1981.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/0377-0257(81)80020-1.
  * Dekée et al. (1986) D. Dekée, P.J. Carreau, and J. Mordarski.  Bubble velocity and coalescence in viscoelastic liquids.  _Chemical Engineering Science_ , 41(9):2273 – 2283, 1986.  ISSN 0009-2509.  doi: http://dx.doi.org/10.1016/0009-2509(86)85078-3.
  * Gummalam et al. (1988) S. Gummalam, K.A. Narayan, and R.P. Chhabra.  Rise velocity of a swarm of spherical bubbles through a non-newtonian fluid: Effect of zero shear viscosity.  _International Journal of Multiphase Flow_ , 14(3):361 – 373, 1988.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/0301-9322(88)90050-X.
  * Kee et al. (1990) D. De Kee, R.P. Chhabra, and A. Dajan.  Motion and coalescense of gas bubbles in non-newtonian polymer solutions.  _Journal of Non-Newtonian Fluid Mechanics_ , 37(1):1 – 18, 1990.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/0377-0257(90)80001-G.
  * Manjunath and Chhabra (1992) M. Manjunath and R.P. Chhabra.  Free rise velocity of a swarm of spherical gas bubbles through a quiescent power law liquid.  _International Journal of Engineering Science_ , 30(7):871 – 878, 1992.  ISSN 0020-7225.  doi: http://dx.doi.org/10.1016/0020-7225(92)90016-A.
  * Miyahara and Yamanaka (1993) Toshiro Miyahara and Shuichi Yamanaka.  Mechanics of motion and deformation of a single bubble rising through quiescent highly viscous newtonian and non-newtonian media.  _Journal of Chemical Engineering of Japan_ , 26(3):297–302, 1993.  doi: 10.1252/jcej.26.297.
  * Liu et al. (1995) Y. J. Liu, T. Y. Liao, and D. D. Joseph.  A two-dimensional cusp at the trailing edge of an air bubble rising in a viscoelastic liquid.  _Journal of Fluid Mechanics_ , 304:321–342, 12 1995.  ISSN 1469-7645.  doi: 10.1017/S0022112095004447.  URL http://journals.cambridge.org/article_S0022112095004447.
  * Rodrigue et al. (1996) D. Rodrigue, D. De Kee, and C.F. Chan Man Fong.  An experimental study of the effect of surfactants on the free rise velocity of gas bubbles.  _Journal of Non-Newtonian Fluid Mechanics_ , 66(2):213 – 232, 1996.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/S0377-0257(96)01486-3.
  * Li et al. (1997) H.Z. Li, Y. Mouline, L. Choplin, and N. Midoux.  Chaotic bubble coalescence in non-newtonian fluids.  _International Journal of Multiphase Flow_ , 23(4):713 – 723, 1997.  ISSN 0301-9322.  doi: http://dx.doi.org/10.1016/S0301-9322(97)00004-9.
  * Chhabra (1998) R. P. Chhabra.  Rising velocity of a swarm of spherical bubbles in power law fluids at high reynolds numbers.  _The Canadian Journal of Chemical Engineering_ , 76(1):137–140, 1998.  ISSN 1939-019X.  doi: 10.1002/cjce.5450760118.  URL http://dx.doi.org/10.1002/cjce.5450760118.
  * Gauri and Koelling (1998) Vishal Gauri and W. Kurt Koelling.  The motion of long bubbles through viscoelastic fluids in capillary tubes.  _Rheologica Acta_ , 38(5):458–470, 1998.  ISSN 1435-1528.  doi: 10.1007/s003970050197.  URL http://dx.doi.org/10.1007/s003970050197.
  * Dubash and Frigaard (2007) N. Dubash and I.A. Frigaard.  Propagation and stopping of air bubbles in carbopol solutions.  _Journal of Non-Newtonian Fluid Mechanics_ , 142(1):123 – 134, 2007.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/j.jnnfm.2006.06.006.  URL http://www.sciencedirect.com/science/article/pii/S0377025706001637.
  * Tsamopoulos et al. (2008) J. Tsamopoulos, Y. Dimakopoulos, N. Chatzidai, G. Karapetsas, and M. Pavlidis.  Steady bubble rise and deformation in Newtonian and viscoplastic fluids and conditions for bubble entrapment.  _Journal of Fluid Mechanics_ , 601:123–164, 2008.  doi: 10.1017/S0022112008000517.  URL http://www.scopus.com/inward/record.url?eid=2-s2.0-42449160607&partnerID=40&md5=c188f7b8fcb9c9f30299e32bba9d8c77.
  * Sikorski et al. (2009) Darek Sikorski, Herve Tabuteau, and John R. de Bruyn.  Motion and shape of bubbles rising through a yield-stress fluid.  _Journal of Non-Newtonian Fluid Mechanics_ , 159(1):10 – 16, 2009.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/j.jnnfm.2008.11.011.  URL http://www.sciencedirect.com/science/article/pii/S0377025708002140.
  * Kishore et al. (2013) Nanda Kishore, V. S. Nalajala, and Raj P. Chhabra.  Effects of contamination and shear-thinning fluid viscosity on drag behavior of spherical bubbles.  _Industrial & Engineering Chemistry Research_, 52(17):6049–6056, 2013.  doi: 10.1021/ie4003188.
  * Tripathi et al. (2015) Manoj Kumar Tripathi, Kirti Chandra Sahu, George Karapetsas, and Omar K. Matar.  Bubble rise dynamics in a viscoplastic material.  _Journal of Non-Newtonian Fluid Mechanics_ , (0):article in press, 2015.  ISSN 0377-0257.
  * Kulkarni and Joshi (2005) Amol A. Kulkarni and Jyeshtharaj B. Joshi.  Bubble formation and bubble rise velocity in gas-liquid systems: A review.  _Industrial & Engineering Chemistry Research_, 44(16):5873–5931, 2005.
  * Chhabra (2006) R.P. Chhabra.  _Bubbles, Drops, and Particles in Non-Newtonian Fluids_.  CRC Press, 2006.
  * Abbassi et al. (2018) W. Abbassi, S. Besbes, M. Elhajem, H. Ben Aissia, and J.Y. Champagne.  Numerical simulation of free ascension and coaxial coalescence of air bubbles using the volume of fluid method (vof).  _Computers & Fluids_, 161:47–59, 2018.  ISSN 0045-7930.
  * Gumulya et al. (2017) M. Gumulya, R.P. Utikar, G.M. Evans, J.B. Joshi, and V. Pareek.  Interaction of bubbles rising inline in quiescent liquid.  _Chemical Engineering Science_ , 166:1–10, 2017.  ISSN 0009-2509.
  * Watanabe and Sanada (2006) Masao Watanabe and Toshiyuki Sanada.  In-line motion of a pair of bubbles in a viscous liquid.  _Jsme International Journal Series B-fluids and Thermal Engineering_ , 49:410–418, 2006.
  * Yuan and Prosperetti (1994) H. Yuan and A. Prosperetti.  On the in-line motion of two spherical bubbles in a viscous fluid.  _Journal of Fluid Mechanics_ , 278:325–349, 1994.
  * Legendre et al. (2003) Dominique Legendre, Jacques Magnaudet, and Guillaume Mougin.  Hydrodynamic interactions between two spherical bubbles rising side by side in a viscous liquid.  _Journal of Fluid Mechanics_ , 497:133–166, 2003.
  * Yu et al. (2011) Zhao Yu, Hui Yang, and Liang-Shih Fan.  Numerical simulation of bubble interactions using an adaptive lattice boltzmann method.  _Chemical Engineering Science_ , 66(14):3441–3451, 2011.  ISSN 0009-2509.  10th International Conference on Gas–Liquid and Gas–Liquid–Solid Reactor Engineering.
  * Katz and Meneveau (1996) J. Katz and C. Meneveau.  Wake-induced relative motion of bubbles rising in line.  _International Journal of Multiphase Flow_ , 22(2):239–258, 1996.  ISSN 0301-9322.
  * Ruzicka (2000) M.C Ruzicka.  On bubbles rising in line.  _International Journal of Multiphase Flow_ , 26(7):1141–1181, 2000.  ISSN 0301-9322.
  * 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.  ISSN 00219991.
  * 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.
  * 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 et al. (2015b) 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, 2015b.
  * 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.
  * 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.
  * Vanka et al. (2016) S. P. Vanka, P. Kumar, K. Jin, and B. G. Thomas.  Single and multiphase flow computations on graphics processing units.  _International Conference on Computational Methods for Thermal Problems_ , (217349), 2016.  ISSN 2305-5995.
  * 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. (2012) Jeremy Horwitz, Purushotam Kumar, and Pratap Vanka.  Simulations of multiphase flow in a t-junction and distributor header.  _Bulletin of the American Physical Society_ , 57, 2012.
  * Horwitz et al. (2013) Jeremy Horwitz, Purushotam Kumar, and Pratap Vanka.  Simulations of three-dimensional droplet deformation in a square-duct at moderate reynolds numbers.  _Bulletin of the American Physical Society_ , 58, 2013.
  * Kumar et al. (2013) Purushotam Kumar, Jeremy Horwitz, and Surya Vanka.  A three-dimensional numerical study of immiscible droplet deformation in a right angle bend.  _Bulletin of the American Physical Society_ , 58, 2013.
  * 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.
  * Horwitz et al. (2014) JAK Horwitz, P Kumar, and SP 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\.
  * 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.
  * Vanka et al. (2015) Surya Pratap Vanka, Kai Jin, Purushotam Kumar, and Brian Thomas.  Rise of an argon bubble in liquid steel in the presence of a transverse magnetic field.  _Bulletin of the American Physical Society_ , 60, 2015.
  * Kumar and Vanka (2022) Purushotam Kumar and Pratap Vanka.  Dynamics of Argon Gas Bubble Pair Rising in Liquid Steel in the Presence of a Transverse Magnetic Field.  _Bulletin of the American Physical Society_ , 67, 2022.
  * Kumar and Vanka (2023) Purushotam Kumar and Surya Pratap Vanka.  Dynamics of argon gas bubbles rising in liquid steel in the presence of transverse magnetic field.  _arXiv preprint arXiv:2301.01797_ , 2023.
  * Vanka et al. (2016) Pratap Vanka, Purushotam Kumar, and Kai Jin.  Numerical Simulation of Turbulent Bubbly Flow in a Vertical Square Duct.  _Bulletin of the American Physical Society_ , 61, 2016.
  * Kumar and Vanka (2021) Purushotam Kumar and Surya Pratap Vanka.  Large scale gpu based simulations of turbulent bubbly flow in a square duct.  _arXiv preprint arXiv:2104.01636_ , 2021.
  * Ohta et al. (2003) Mitsuhiro Ohta, Eiji Iwasaki, Eiji Obata, and Yutaka Yoshida.  A numerical study of the motion of a spherical drop rising in shear-thinning fluid systems.  _Journal of Non-Newtonian Fluid Mechanics_ , 116(1):95 – 111, 2003.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/j.jnnfm.2003.08.004.
  * Zhang et al. (2010) Li Zhang, Chao Yang, and Zai-Sha Mao.  Numerical simulation of a bubble rising in shear-thinning fluids.  _Journal of Non-Newtonian Fluid Mechanics_ , 165(11-12):555 – 567, 2010.  ISSN 0377-0257.  doi: http://dx.doi.org/10.1016/j.jnnfm.2010.02.012.
  * Ohta et al. (2012) Mitsuhiro Ohta, Sachika Kimura, Tomohiro Furukawa, Yutaka Yoshida, and Mark Sussman.  Numerical simulations of a bubble rising through a shear-thickening fluid.  _Journal of Chemical Engineering of Japan_ , 45(9):713–720, 2012.  doi: 10.1252/jcej.12we041.
  * Carreau (1972) Pierre J. Carreau.  Rheological equations from molecular network theories.  _Transactions of The Society of Rheology_ , 16(1), 1972\.