The Influence of Bubbles on Foamed Cement Viscosity Using an Extended Stokesian Dynamics Approach

Abstract

We want to study the influence of bubbles on the viscosity of suspensions with a computational approach that also accounts for the arrangement of the bubbles due to shearing flow. This requires a large number of bubbles to properly simulate and requires a large amount of computational resources. Here we develop a set of equations to define the viscosity ratio from the simulation results to show the influence of the bubbles on the viscosity as a function of the volume fraction. One application of this work has been used to study a specific type of cement that has bubbles injected into the slurry while it is still fluid. The bubbles are added to reduce the density but they also improve the properties of the cement with the increase in viscosity. We show that the computed results match the few experimental results that have been reported.

1. Introduction

The primary motivation for this work is to study a specialized wellbore cement that has bubbles dispersed and suspended in it. The cement is referred to as foamed cement and the cement is foamed to reduce its density while still maintaining compressive strength. The cement density is designed to keep the cement gradient between the fracture and pore pressures [1,2] and current technologies allow the foamed cement density or quality (gas volume fraction) to be controlled and optimized during operations [3]. Foamed cements offer an economically viable low-density option that still maintains compressive strength. Foaming the cement also offers other improved properties compared to conventional cement or alternative light-weight cement [1,4,5,6,7].

In drilling and completing a well, cement is placed in the annulus between the well casing and the formation to prevent fluid migration between the formation layers and to prevent fluids from migrating to the surface. The cement also provides support to the steel casing pipe in the center of the well. A successful cement job will provide complete “zonal isolation” [8,9,10], where fluids in the wellbore remain isolated. Any means for fluids to migrate through the cement via fluid pathways is undesirable and could lead to severe issues.

Foamed cement is typically used in wells with weak formations or anywhere a low density cement is required but it has been used successfully in other adverse environments where the foaming provides better performance, over conventional or unfoamed cements [11]. It has been used for high-temperature/high-pressure wells where stress cracking is of concern due to the cycling pressures and temperatures [4,12]. Its higher ductility when compared to conventional cements makes it more resilient to stress cracking [1,4,5,7]. With its ability to withstand cyclical conditions, foamed cement is well suited for geothermal wells [5]. Foamed cement reduces lost circulation [5,8], where some of the cement flows into the formation, since foamed cement can be maintained below the formation fracture gradient during placement. Foamed cement offers improved mud displacement and controls gas migration better than conventional cements [1,5,7] due to the presence of the bubbles in the cement slurry. Even with the specialized foaming equipment, overall costs can be reduced by using foamed cement and, due in part to using less material, foamed cement has a lower environmental impact [5].

Foamed cement is ‘designed’ for the downhole conditions where it is to be placed in the well [11]. The foam is created when the base cement slurry flows into a tee met at a right angle by pressurized nitrogen, which is dispersed into the cement slurry via some form of atomizer. Stabilizers and additives are generally pumped into the cement slurry prior to foaming [4]. The design of the cement is based on the midpoint of cement placement depth so the changing hydrostatic pressure gradient and temperature changes occurring as the foamed cement is placed in the well must also be considered in the final cement design [3].

The American Petroleum Institute (API) recommends a foam quality (the volume fraction of added gas) below 35% for a stable foam that will maintain mechanical integrity and proper zonal isolation [3]. Above 35% volume of gas in the cement, the permeability may increase to unacceptable levels, and the compressive strength and ability to support the casing could be compromised. The foam should also have well dispersed bubbles that remain suspended during placement in the well without significant clustering, coalescence, or other undesirable configuration changes until the cement fully hardens [11]. It is thought that the bubble sizes and extent to which the bubbles remain well dispersed affects the overall strength and performance of the cement. Coalescence or clustering of the bubbles could provide migration pathways for fluids once the cement is fully hydrated and could reduce the compressive strength of the cement sheath.

2. Methods

2.1. Model Assumptions

We used a computation method to calculate the viscosity ratio (also called the relative viscosity to describe the change in viscosity due to the suspended bubbles). The computation method that we developed [13] extends Stokesian Dynamics and Fast Lubrication Dynamics. To employ the computation methods, several assumptions must be made. The bubbles are assumed to be spherical and discrete and remain spherical throughout the simulation. Foamed cement samples that have been produced in the laboratory and collected as part of a collaborative field effort show spherical, stable, discrete bubbles [14,15]. Though the interaction considers spherical bubbles, the surface is allowed to deflect in the region of bubble pair interaction [16]. Another assumption is that the bubbles are neutrally buoyant. This has also been observed experimentally [14,15,17] and were observed to remain in place during curing under pressure [15]. We also must assume that the suspending fluid is Newtonian. Cement is known to have a yield stress and nonlinear behavior but this assumption must be made to maintain linear equations.

2.2. Stokesian Dynamics and Fast Lubrication Dynamics

Stokesian Dynamics [18] and Fast Lubrication Dynamics (FLD) [19,20] include both close-range and long-range interactions to simulate spherical objects suspended in a base fluid. They both employ methods to incorporate these ranges of interactions in more efficient dynamics so that large numbers of spherical objects can be simulated at any volume fraction of objects. Even in dilute suspensions, long-range effects need to be considered, because in a suspension, the motion of each suspended object (to also include bubbles) is transmitted through the base fluid. In the quasi-static (creeping flow) limit, this is felt immediately throughout the entire system by all the other objects. For a dilute suspension of objects, the objects are far apart so the detailed shape and structure of the suspended objects does not matter to leading-order. The velocity disturbance of a suspended object decays like a point force as $1/r^2$, where r is the radial distance from the object [18]. When the bubbles are close together, the interaction of each pair is dominated by a pairwise force which comes from lubrication theory, with specific surface properties representative of a bubble [21]. The Stokesian Dynamics method accounts for the far-field interaction through multipole expansions and for the near-field interactions through pairwise interactions [18]. The FLD method further increases the efficiency by using fast approximate methods for the far-field interaction [19,20]. We employ an extended FLD method with specific pairwise near-field interactions for bubbles [16].

2.3. Governing Equations

For the Stokes flow regime, the viscous forces dominate and the inertial forces have negligible effects. The movement of a particle is governed by the equation [22,23]:

(1)$\mathit{m}·\frac{\mathrm{d}\mathit{U}}{\mathrm{d}t}={\mathit{F}}^H$

Here, $\mathit{U}$ is the generalized velocity vector of the particles with entries corresponding to both linear and angular velocities. The velocity vector therefore has 6 components, corresponding to linear and angular velocities along the Cartesian directions, for each particle. Correspondingly, $\mathit{m}$ is the generalized moment of inertia matrix containing entries for mass as well as the rotational inertia and, $\mathit{F}$ is the generalized force vector with entries for both force as well as moment. Contained in the hydrodynamic force vector, ${\mathit{F}}^H$, are then the 6 Cartesian components of the force and torque for each particle.

The hydrodynamic forces can be determined from the generalized equation relating the velocities and particle positions, $\mathit{X}$, to the forces.

(2)$\left(\begin{array}{c}{\mathit{F}}^H\end{array}\right)=\mathcal{R}\left(\mathit{X}\right)·\left(\begin{array}{c}\mathit{U}\end{array}\right)$

$\mathcal{R}$ is the linear resistance matrix that is a function of the particle positions, $\mathit{X}$. $\mathcal{R}$ is almost the inverse of the mobility matrix, $\mathcal{M}$, and is a close approximation, though they are not complete inverses of each other [21]. The mobility matrix is more computationally efficient to construct but does not capture the physics of the problem as well as the resistance matrix and can allow the spherical objects to overlap in an unrealistic way. Stokesian Dynamics takes advantage of these properties to reduce computation time. Stokesian Dynamics uses the inverse of the grand mobility matrix, ${\left({\mathcal{M}}^{\infty }\right)}^{-1}$, as an approximation to the grand resistance matrix, ${\mathcal{R}}^{\infty }$, which captures the long-range effects. In FLD, the grand mobility matrix is replaced with an isotropic resistance tensor, ${\mathcal{R}}_{Iso}$, to further improve computational efficiency [19,20]. Stokesian dynamics and FLD both explicitly include the close-range two-body lubrication interactions, ${\mathcal{R}}_{2B}$, using the lubrication equations developed by Kim and Karrila [21], Jeffrey and Onishi [24] and given succinctly and used in the form developed in Ball and Melrose [25].

The resistance matrix used in FLD is described by:

(3)${\mathcal{R}}_{FLD}={\mathcal{R}}_{Iso}+{\mathcal{R}}_{2B}-{\mathcal{R}}_{2B}^{\infty }$

Since the construction of the mobility matrix and inversion (in Stokesian Dynamics) or the construction of the isotropic resistance matrix approximation include all the pair interactions, the two-body close-range interactions, ${\mathcal{R}}_{2B}^{\infty }$, need to be subtracted to not include them twice. Though, part of the ${\left(\mathcal{R}\right)}_{FLD}$ is an approximation, the dominant forces are the lubrication pairwise forces. Ball and Melrose [25] showed that accurate simulations of spherical particles can be performed with only the close-range lubrication interactions and achieve results comparable to more accurate methods for volume fractions of particles of $0.20$ and above.

2.4. Bubble Interaction Modeling

Our approach in simulating bubbles is to consider three properties of bubbles: (1) slip surface conditions (${\mathit{F}}^{Bubble}$); (2) the ability of the surface to deflect locally in the interacting region (${\mathit{F}}^{Elastic}$); and (3) the repulsive and attractive quality of the surface due to a surfactant (${\mathit{F}}^{Direct}$). The total of these three forces serve to include the predominant properties of realistic bubbles. A direct force that represents the surfactant also serves to prevent nonphysical overlaps of the bubbles during the simulations. A depiction of the three forces is shown in Figure 1. The total hydrodynamic force for the bubbles is then defined by [16]:

(4)${\mathit{F}}^{Total}={\mathit{F}}^{Bubble}+{\mathit{F}}^{Elastic}+{\mathit{F}}^{Direct}$

Real bubbles have a mobile interface. Theoretically, the slip on the bubble surface can be considered very high with no shearing resistance. If there is a high amount of slip on the bubble surface then there is no shear traction at the surface boundary condition. Under these conditions, in a pair interaction where the bubbles interact with a relative shearing motion, the bubbles will just slip by each other [21,26]. It can then be assumed that the bubbles have no shearing resistance so the shearing resistance terms all become zero. To simulate surfactant stabilized bubbles though, we need to have a mobile surface boundary and an additional realistic normal direct force to prevent bubble overlap as the surfactant does in real bubble suspensions. The force on approaching bubbles is weaker and will allow bubbles to touch in finite time so we must impose a direct force to avoid unreal overlaps. The direct force was chosen so that it enabled the bubble pairs to approach very close to each other to allow a small surface deflection but did not allow complete bubble overlap.

The hydrodynamic interaction for a pair of bubbles, with the same size diameter, has been derived by [21] for the case of two bubbles approaching each other and is given by [16,21]:

(5)${\mathit{F}}^{Bubble}=a_{sq}^{Bubble}\mathit{U}$

where $a_{sq}^{Bubble}$ is [21]:

(6)$a_{sq}^{Bubble}=\pi {\mu }_{f}a\left(2log\frac{a}{h}+4\left({\gamma }_{Euler}+log2\right)\right)$

${\mu }_f$ is the viscosity of the suspending fluid, a is the radius of the bubble, and h is the gap between the bubbles. ${\gamma }_{Euler}$ is Euler’s constant (${\gamma }_{Euler}$ = 0.577216…).

The elastic force used in the bubble interaction is taken directly from [27]:

(7)${\mathit{F}}^{Elastic}=a_{Elas}E\delta$

(8)$a_{Elas}=\sqrt{a(2h+\left|\delta \right|)}$

(9)$\delta =\frac{(2a+h-r)}{2}$

In the region of interaction between the bubble pair, the surface is considered to become flattened. $a_{Elas}$ is the region of the bubble surface that is interacting in the lubrication force after deflection has caused the surface to become flattened. $\delta$ is the deflection of the bubble surface. Once the bubbles move away from each other, the deflection of the local surface region is zero and the spherical shape of the bubble is restored. To limit the amount that the bubble surface can deflect, a maximum deflection value was set. E represents an effective modulus to deform the interface and bubble and was set to ${10}^2$ for the results shown here. The value of E could be adjusted to match experimental data.

The direct force that was chosen and used for the simulation results reported here is the Lennard-Jones potential force. Surfactants can be attractive and repulsive in nature [28], which is why the Lennard-Jones force was chosen as the direct force. The surfactant forces may include van der Waals and hydrostatic attraction, electrostatic repulsion, steric repulsion, etc. [28], however, it is generally not possible to obtain the specific details about the surfactant due to its proprietary properties. So, no specific surfactant was considered. Using a generalized repulsive force that increases steeply when the particles come in close relative positions to each other has been done in Soft Dynamics [27,29] and the direct force used here is deployed in a similar manor. The Lennard-Jones potential included with LAMMPS [30] was used as the direct force. The Lennard-Jones force in LAMMPS is:

(10)${\mathit{F}}^{Direct}=\frac{4ϵ}{\sigma }\left[9{\left(\frac{\sigma }{r}\right)}^{10}-6{\left(\frac{\sigma }{r}\right)}^7\right]$

This represents a Lennard-Jones potential that is less repulsive than the standard Lennard-Jones potential (The standard Lennard-Jones potential is ${\varphi }_{LJ}=4ϵ\left[{\left(\frac{\sigma }{r}\right)}^{12}-{\left(\frac{\sigma }{r}\right)}^6\right]$). $ϵ$ is the energy scale and $\sigma$ is the length scale. The influence of the values of $ϵ$ and $\sigma$ on system dynamics are described in detail in [16], so we do not present them again here. A value of $0.008$ was used for $ϵ$ and a value of $0.9$ was used for $\sigma$. ${\left(\frac{\sigma }{r}\right)}^{10}$ is the repulsive part and ${\left(\frac{\sigma }{r}\right)}^7$ is the attractive portion.

To simulate bubbles of different sizes, a reduced radius [21,31] was used. If the radius of particle a is a and the radius of particle b is b, then the reduced radius is:

(11)$a_{reduced}=\frac{a\times b}{a+b}$

To simulate polydisperse bubbles, $a_{reduced}$ is used in the equations instead of a.

2.5. Simulation Framework

LAMMPS (lammps.sandia.gov Open-source classical (non-quantum) molecular dynamics code developed and maintained at the Sandia National Labs.) (Large-scale Atomic/Molecular Massively Parallel Simulator) was modified with the interactions described in this paper and was used to run the simulations.

2.6. Simulation Inputs

Simulation inputs (i.e., the bubble size and each bubble’s $x-,y-,z-$position) were created by randomly placing bubbles with a diameter of 1 in a simulation box of size $10\times 10\times 10$. Volume fractions of $10\%$, $20\%$, $30\%$, $40\%$, $45\%$, and $50\%$ were created with increasing numbers of bubbles in the simulation box. To remove any initial overlapping bubbles, a soft potential was used to push bubbles apart and remove overlaps. Through a series of test simulations, it was determined that replicating the simulation box four times in each direction eliminated system size effects on the relative viscosity. However, due to replicating the same simulation box, the bubble positions were repeated in each box. To add back the random placement of the bubbles, a Brownian pair interaction was used to randomly move the bubbles a small amount and re-introduce the random nature back into the initial bubble positions. These steps were not part of the dynamics of the simulations but served to create the initial inputs to the simulations.

To create polydisperse bubble systems, the monodisperse bubble positions were used and the bubble size was increased and decreased randomly on the monodisperse bubbles to add slight polydispersity to the distribution of bubble sizes. The system of polydisperse bubbles were then replicated and the positions randomized as described for monodisperse bubbles.

2.7. Implementation of the Shearing Flow

The shearing flow was implemented as described in [13,16,32]. Lees-Edwards boundary conditions were imposed (See the LAMMPS User Manual [30] for the implementation of the Lees-Edwards like boundary conditions, which uses a function called “fix deform” to apply a strain rate to the simulation box in the specified orthogonal box directions.). The strain was imposed on the simulation box in the $xy$-direction so that the velocity of each bubble becomes a function of its position in the y-direction. The periodic boundary conditions were imposed such that when a bubble crossed the simulation box boundary in any direction, the velocity of the bubble was remapped to correspond to the new position in the simulation box (See e.g., Dayal and James [33] for a discussion of this.). The simulations were all sheared until the stress reached a constant value and the value of $\dot{\gamma }\times t_{total}$ reached 200 so that the simulation results could all be compared. The length of time that the simulation was run is $t_{total}$ and $\dot{\gamma }$ is the strain rate.

An explicit time-integration scheme was used for these interactions and therefore required a small time step. Testing was performed to optimize the dimensionless time step, $\Delta t$. It was determined that the time step should be kept below a value of $0.002$ for convergence. For the results contained herein, a time step of $\Delta t=0.001$ was used. Smaller time steps would require the simulation to run for a higher total time to reach the same value of total strain.

2.8. Calculation of Viscosity Ratio of Simulated Suspensions

The relative viscosity can be calculated from the summation of the stress on each bubble. The stress on each bubble is a function of the force, $\mathit{F}$, on each bubble due to its $N_p$ neighbors. The stress is calculated in LAMMPS at each time step and is given by [30,34]:

(12)${\sigma }_{i,\alpha \theta }=\frac{1}{V_i}[-\sum _{n=1}^{N_p}\left(r_{1_{\alpha }}{F}_{1_{\theta }}+r_{2_{\alpha }}{F}_{2_{\theta }}\right)]$

$\alpha$ and $\theta$ run over the coordinate directions to compute the 6 components of the symmetric stress tensor. ${\mathit{r}}_1$ and ${\mathit{r}}_2$ are the positions of bubble 1 and bubble 2, respectively in the pair that has pairwise interactions. And, ${\mathit{V}}_{\mathit{i}}$ is the volume of each bubble. In the typical calculation, the value does not consider the bubble volume so to properly calculate the relative viscosity, the total volume, $\mathit{V}$, of the simulation box must be included in the relative viscosity calculation. The relative viscosity or stated another way, the viscosity ratio (i.e., $\mu$ over the viscosity of the suspending fluid, ${\mu }_f$), can be calculated from the average of the total stress once the system has reached steady state:

(13)${\mu }_{\mathrm{relative}}=\frac{\mu }{{\mu }_f}=\frac{{\sum }_i{\sigma }_{i,xy}}{\dot{\gamma }{\mu }_{f}V}$

3. Results

3.1. Foamed Cement Viscosity from Prior Experimental Results

3.1.1. Cement Rheology

Cement slurry rheology is an important parameter to understand for designing and implementing a cement job [35,36]. Cement slurries are a suspension of small particles undergoing chemical reactions and are known to have Non-Newtonian rheology, generally assumed to follow a Bingham model [35,37]. Bingham fluids have a yield stress below which no flow occurs and after the yield stress is attained, the fluid maintains a constant viscosity. Cement slurry rheology, however, can be time dependent and shear dependent and depend on shearing history and resting time [35,37]. Other models that have been used to describe cement slurry rheology are the power law model, with limitations, and the Herschel-Bulkley model [35]. Both of these models account for the shear thinning behavior of cement and the Herschel-Bulkley model includes a yield stress. The time dependent behavior of rheology—thixotropy: time-dependent shear thinning or rheoplexy: time-dependent shear thickening—is considered to be reversible. With cement, however, the hydration process will also change the rheology as time progresses [35,37] so it is not a truly reversible process. Cement slurry rheology is also dependent on shearing and resting history. Particles within the suspension can cluster together due to random movements under static or no flow conditions [35,37,38]. During shearing, the particle clusters can be broken apart and this may be the mechanism of shear thinning that occurs with time [35,37,38]. At low shear rates, cement slurries have been shown to be thixotropic, which means that as they are sheared, the viscosity decreases with time, though thickening can occur due to the hydration reaction. Foamed cement is a complex system and inherently dynamic, which makes it a challenge to study rheology with reliable results. Cement rheology is a function of several factors including the cement slurry design and the cement slurry production methods, and can be affected by the measurement techniques themselves [35].

3.1.2. Rheology of Foamed Cement

Foamed cement rheology is further complicated by the addition of bubbles suspended in the cement. The changing downhole conditions in the well as the cement is pumped can affect the rheology as the quality changes. Few studies have actually been done to explore the rheological properties of foamed cement, even at ambient conditions. In the book Well Cementing [10], Chapter 14: Foamed Cement [39] a comment is made about the lack of rheological information and the importance of the data. At the time only one study had been referenced [40] and only two other studies [41,42] have been done since the publishing of that book. Some of the conclusions of the studies contradict each other. The earliest study [40] uses two different capillary flow tubes (0.091 inch I.D. × 8 ft long, 0.130 inch I.D. × 8 ft long) to measure the rheology. A description of the cement, even the type used, was not included in this thesis. The cement mixer is not described fully and these specifics can influence the rheology. The cement was mixed with a “cement mixer” and injected into the system with a “variable flow rate, positive displacement plunger pump”, rated to 3000 psi. The foaming agent was injected into the base slurry with a “variable flow rate positive displacement gear pump” and mixed into the cement as the slurry flowed through a small tube coil. After air injection, the cement slurry was foamed by passing it through three stainless steel foam generator tubes (each 12 inch long × 0.33 inch I.D.) that were filled with stainless steel shavings to cause mixing and foaming.

Results span foamed cement qualities from 30% up to 65% and shear rate values from 1000 1/s to 10,000 1/s [40]. It was shown that the viscosity increased with increasing quality for a given shear rate and the viscosity decreased with increasing shear rate. At high shear rates ( 7000 1/s and higher) and qualities lower than 50%, the viscosity seems to become constant, however, typical shear rates during placement are below 1000 1/s [39]. The viscosity as a function of quality varied between the two diameter flow tubes, with the larger tube having higher reported values. They conclude that the Power Law model or the Bingham plastic model describe the data. While a first of its kind study, pertinent details were left out of this thesis and the shear rates explored are well above those encountered during a cementing job.

A later study [41] uses an inline flow-through rotational viscometer to measure the rheology. The base slurry consisted of Class G cement, water, and a few additives. The base slurry was made by first mixing the water, dispersant, and a fluid loss additive at 200 rpm for about 1 min. Then the cement was added while mixing at 200 rpm for 50 s. The cement slurry (without foamer) was agitated at 1500 rpm for 5 min. Rheology measurements were taken using a regular rotational bob and sleeve viscometer. The original base slurry was removed from the viscometer, placed back in the original mixing container, foamer was added and mixed gently without foaming it. To foam the cement, the cement mixture with foamer was injected into a mixing chamber where pressurized nitrogen was also injected. The cement, additives, and foamer were agitated at 1600 rpm for 5 min with the mixer. To measure viscosity, the foamed cement flowed through the inline viscometer at a constant rate and measured at various shear rates.

The study of [41] is the first and only known attempt to assess the influence of pressure on the rheology. They also are the first to report viscosity values within relevant shear rates (1000 1/s or less). In conclusion, they find that a Power Law model provides the best fit to the data but at lower shear rates. However, the cement slurry has a yield stress, which is not captured with a Power Law model. They also show that the viscosity increases with quality for 10%, 20%, and 30% foam qualities. The discrepancy with previous reports/experience comes in their conclusion that the base slurry viscosity is higher than both the 10% and 20% foam qualities measured and, rather, is close to that of the 30% foam quality and is almost the same as the 30% foam quality for shear rates of 50 1/s and above. They question the claims that foamed cement displaces drilling muds better than neat cement as experience from industry suggests. It is unclear if these results are real or due to experimental issues. It is also possible that in measuring the foam viscosity, the foam was destabilized and the foamer reduced the slurry viscosity as compared to the ‘base’ slurry, which had no foamer. The previous study from 1977 does not report data for qualities below 30% or for the neat cement.

The most recent study, [42], uses a uniquely shaped paddle rather than a bob and sleeve viscometer to measure the viscosity and yield point. The paddle devise, a Fann Yield Stress Adapter (FYSA), was used to minimize wall slippage that can be an issue with bob and sleeve viscometers [42]. The FYSA is designed to minimize wall slippage, measures higher volumes of material, and keeps the bubbles and cement particles better dispersed. The base slurry was prepared with Class H cement and water. After initial mixing at high speed in a blender, the foaming agent was added and stirred by hand without producing any foaming and then mixed in a blender to foam as recommended according to API Recommended Practices [43]. The viscosity measurements were then conducted in the viscometer unlike the former recorded measurements, which were done with flow through methods. The results can be considered the most accurate viscosity measurements of foamed cement due to the advantages and reduction in errors due to the design of the FYSA. Measurements were reported for a range of volume fractions of gas: 10%, 20%, 30%, and 40%.

It is challenging to measure the rheology of a foamed cement as the methods used could affect the rheology [39]. Foamed cements used in well cementing are pressurized and to replicate pressurized foamed cement in the laboratory requires specialized equipment to maintain the foam under pressure and perform the measurements under pressure.

3.2. Viscosity Ratio of Suspensions

An important result of adding bubbles to a fluid is to increase the effective viscosity. In this section, we examine the effect of the bubble volume fraction in the suspending fluid on the viscosity, calculated from the simulated results. We provide equations to describe the viscosity ratio as a function of the bubble volume fraction. The equations represent a valuable contribution to determine the influence of the bubbles on the viscosity compared to the base fluid alone.

The definition of the viscosity ratio for the simulated results is given by Equation (13). This equation was used to compute the viscosity ratio determined from each simulation that was run. The values of the computed viscosity ratio can then be fitted to equations and the coefficients determined to define the relationship between bubble volume fraction and the subsequent increase in viscosity.

The relationship between the volume fraction ($\varphi$) of objects in a Newtonian fluid (with a viscosity of ${\mu }_f$) and the corresponding change in viscosity can be described by several equations with varying levels of complexity. The effects of object concentration on the suspension viscosity is summarized in [44]. As the bubbles, in the foamed cement, act like particles, we examine typical equations used to describe suspensions.

Einstein’s classical equation accounts for the increase in viscosity due to the object volume fraction and is linear in form. The Einstein equation is therefore valid for low object volume fractions. The form of the Einstein Equation [45] is:

(14)$\mu ={\mu }_f(k+l\varphi )$

Krieger and Dougherty account for the limits of the amount of objects that can be in the suspension with the maximum packing fraction (${\varphi }_{max}$). The form of the Krieger and Dougherty Equation [46] is:

(15)$\mu ={\mu }_f{\left(m-\frac{\varphi }{{\varphi }_{max}}\right)}^{-n{\varphi }_{max}}$

Batchelor and Green’s equation is a second order polynomial and fits the simulation and experimental data the best of all the model equations. The form of the Batchelor and Green Equation [47] is:

(16)$\mu ={\mu }_f(o+p\varphi +q{\varphi }^2)$

The coefficients for the above equations for both particles and bubbles are summarized in Table 1.

Taylor’s equation allows for fluid spheres of a different fluid from the base suspending fluid. However, the Taylor Equation did not fit the simulation results well and so the coefficients are not given. Here, for bubbles, the fluid is air and the viscosity (${\mu }_{bubble}$) is much lower than the suspending fluid. The form of Taylor Equation [48] is:

(17)$\mu ={\mu }_f\left[kk+ll\varphi \left(\frac{{\mu }_{bubble}+mm{\mu }_f}{{\mu }_{bubble}+{\mu }_f}\right)\right]$

The fit of Equations (14)–(16) to the simulated data results are shown in Figure 2 and Figure 3. It is obvious that the second order polynomials, including the Batchelor form of the viscosity equation, fit the data the best. Figure 2 shows the equation fits to the simulated data of the monodisperse spherical rigid particles and monodisperse bubbles. Figure 3 shows the equation fits to the simulated data of the polydisperse spherical rigid particles and polydisperse bubbles. In both figures, the Einstein equation is only fit through the point $[0,1]$ and the first two data points, since it is valid at low volume fractions.

4. Discussion

When comparing the data between the particle suspensions and bubble suspensions, for the case of the same monodisperse particles and bubbles with different surface conditions, the results on the effects to the viscosity are not significant. However, when comparing polydisperse suspensions, which are more realistic, the bubbles do not increase the relative viscosity as much as particles. The viscosity is a function of the stress on each particle or bubble and since the particles are rigid and have no slip surfaces, the overall system stress can increase more at higher volume fractions. The bubbles are able to slide by each other and therefore the overall stress does not increase as much.

The simulated relative viscosity data can be compared with some of the known experimental data (Figure 4). The most comprehensive collection of the viscosity of foamed cement is contained in Al-Mashat’s thesis [40]. However, Al-Mashat’s method for measuring the viscosity is dependent on the size of the pipe used in the setup. So it can be expected that the measurements would be different with another size of pipe. The data does show the increase in the viscosity due to the addition of bubbles. Due to the discrepancy of the results of Ahmed [41] compared to the other experimental results, they are not displayed here along with the simulated results. Olowolagba [42] reports more accurate experimental foamed cement viscosity results and are used to verify the simulated results. The line fitted through the experimental data of Olowolagba [42] is a second order polynomial fit through the data at the lowest shear rate. The simulated data of bubbles suspended in a fluid match the most accurate experimental measurements [42].

5. Conclusions

We have shown that simulated bubble suspensions match with experimental measurements. To determine what effects bubbles will have on a suspension viscosity, the equations with the coefficients determined here can be used to predict the effects of the addition of particles and bubbles to the suspension fluid. Experiment setups can influence or change the properties trying to be measured and require specific and careful consideration and relevant simulations can take significant computing resources. However, the relationships shown herein can be used to predict the influence of bubbles and particles on the viscosity.

Author Contributions

The conceptualization, all authors; methodology, all authors; simulations, E.R.; validation, all authors; formal analysis, E.R.; investigation, all authors; resources, all authors; writing—original draft preparation, all authors; writing—review and editing, all authors; visualization, E.R.; supervision, M.M, K.D.; project administration, E.R.; funding acquisition, E.R., K.D.

Funding

This project was funded by the Department of Energy, National Energy Technology Laboratory an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof; Computing resources were provided by the National Science Foundation through the Extreme Science and Engineering Discovery Environment (XSEDE) program; specifically the XSEDE Bridges Computing Resources at the Pittsburgh Supercomputing Center [grant number TG-DMR160018]; Kaushik Dayal was supported through an appointment to the National Energy Technology Laboratory Faculty Research Participation Program sponsored by the U.S. Department of Energy and administered by the Oak Ridge Institute for Science and Education.

Acknowledgments

Kaushik Dayal acknowledges an appointment to the National Energy Technology Laboratory Faculty Research Participation Program sponsored by the U.S. Department of Energy and administered by the Oak Ridge Institute for Science and Education; Eilis Rosenbaum acknowledges the award of the Bradford and Diane Smith Graduate Fellowship that was created through the generosity of Bradford Smith, a 1969 graduate of the Carnegie Mellon University (CMU) College of Engineering, and his wife Diane, a 1971 CMU graduate of the Dietrich College of Humanities and Social Sciences to help support the graduate studies of highly deserving students in the College of Engineering at CMU; This paper draws from the Ph.D. Thesis of Eilis Rosenbaum at Carnegie Mellon University [32].

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

MDPI	Multidisciplinary Digital Publishing Institute
FLD	Fast Lubrication Dynamics
LAMMPS	Large-scale Atomic/Molecular Massively Parallel Simulator
FYSA	Fann Yield Stress Adapter
API	American Petroleum Institute

Figures and Table

Figure 1. The bubble interaction is a combination of the Bubble Force, FBubble, the Elastic Force, FElastic, and the Direct Force FDirect.

Figure 2. The relative viscosity determined from the simulated results is indicated by points and the lines are the equations fit to the data.

Figure 3. The relative viscosity determined from the simulated results is indicated by points and the lines are the equations fit to the data.

Figure 4. The relative viscosity determined from the simulated results is shown with the experimental data of Al-Mashat [40] and Olowolagba [42].

Table 1

Coefficients for fits of Equations (14)–(17) for the simulated results with a $95\%$ confidence interval in the fit. The equation fit was determined with Matlab.

Coefficient	Monodisperse Particle	Polydisperse Particle	Monodisperse Bubble	Polydisperse Bubble
k	7.463$\times {10}^{-7}$	3.578$\times {10}^{-6}$	2.479$\times {10}^{-8}$	0.4661
l	0.1605	0.1609	0.1665	0.09995
m	1.057	1.014	1.048	0.9535
n	0.02799	0.02725	0.0281	0.0198
o	0.9777	0.9	0.943	1.006
p	0.007433	1.777$\times {10}^{-10}$	3.793$\times {10}^{-14}$	0.01958
q	0.003042	0.003415	0.003393	0.001581