Displacement Flows of Foamed Cement in PrimaryCementing of Oil and Gas wellsbyNikoo Rahimzadeh HanachiB.Sc., Chemical Engineering, Amirkabir University of Technology, 2015A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMaster of Applied ScienceinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Mechanical Engineering)The University of British Columbia(Vancouver)April 2018c© Nikoo Rahimzadeh Hanachi, 2018AbstractThe idea of using foamed cement, which is essentially a bubbly liquid cementslurry, in the procedure of primary cementing of oil and gas wells is to have controlover the density and be able to have lighter cements. However, some complicationsarise from the liquid-gas (compressible) mixture. First, there is an ongoing concernabout the stability of the foamed slurry, e.g. as was queried in the enquiry intothe 2010 BP Macondo incident. Secondly, there are questions regarding how therheology of the foamed slurry should be modelled. Thirdly, there are questionsregarding the stability of the placement flow itself, i.e. assuming that the foamedslurry itself remains intact.Here we have developed two preliminary models of primary cementing thatinclude foamed cements. A one-dimensional hydraulic flow model is derived fordisplacement of foam flow inside the casing and annulus. In the annular regionoutside the casing, we modify the model developed to model laminar displacementflows of incompressible fluid, via a two-dimensional gap-averaged model. Themain differences of our model compared to an incompressible fluid displacementare: (i) flow in our case is represented with a mass-flux stream function. (ii) densityand rheological properties of foam are pressure dependent. Numerical simulationof displacement flows using the resulting model show that the annular flow ex-hibits buoyancy driven instabilities in many situations as the foam advances up theannulus.iiLay SummaryUnderstanding the behavior of cement during its displacement inside the wellboresis very critical in the construction of oil and gas wells. In this work we have devel-oped a model in order to predict the behavior of ” f oamed cement” when it flowsinside a pipe or an annulus. Foamed cement, which is a compressible fluid, con-sists of a base cement slurry (water and dry cement) that suspends a gas phase. Dueto the existence of bubbles, all the properties of the foam are pressure dependent.Hence, as the foam flows down inside the wellbore, due to compression of the gasphase, its density increases. Using the model developed in this work we are ableto predict different properties of the foam at any point inside the wellbore, whichhelps us to figure out how using the foamed cement in oil and gas industry as analternative for the conventional cements changes the efficiency of the displacement.iiiPrefaceThe research presented in this dissertation was conducted by the author, Nikoo Rahimzadeh Hanachi under the guidance of Dr. Ian Frigaard. The MATLAB code that is used in Chapter 5 for solving the two-dimensional mathematical model, was originally developed by I. Frigaard extended by A. Maleki to model turbulent annular displacements. This code has been modified based on the conditions and requirements of this work by the author and with the advice of A. Maleki.Some of the results of this thesis have been submitted for publication as:N. Hanachi, A. Maleki, I.A. Frigaard, A Model for Foamed Cementing of Oil and Gas Wells; submitted February 2018.The work will also be presented at AERC2018, the Annual European Rheology Conference in Sorrento, Italy on April 17-20, 2018: N. Hanachi, A. Maleki, I.A. Frigaard, Displacement of foamed cement in primary cementing, Paper ID: 606.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiNomenclature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 The Oil and Gas Industry in Canada . . . . . . . . . . . . . . . . 11.2 Primary Cementing . . . . . . . . . . . . . . . . . . . . . . . . . 21.3 Why Foamed Cement . . . . . . . . . . . . . . . . . . . . . . . . 51.4 Thesis Objectives and Outline . . . . . . . . . . . . . . . . . . . 62 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1 What Are Foamed Cements . . . . . . . . . . . . . . . . . . . . . 82.2 Rheology of Foamed Cement . . . . . . . . . . . . . . . . . . . . 112.3 Operational Methods in Foamed Cementing Job . . . . . . . . . . 15v2.4 One-Dimensional Hydraulic Models of Foam Flow . . . . . . . . 162.5 Flow in the Annulus . . . . . . . . . . . . . . . . . . . . . . . . . 202.5.1 Two-Dimensional Flow of Incompressible Fluid in Annulus 203 One-Dimensional Version of a General Model . . . . . . . . . . . . . 233.1 One-Dimensional Model . . . . . . . . . . . . . . . . . . . . . . 233.1.1 Wall Shear Stress . . . . . . . . . . . . . . . . . . . . . . 263.1.2 Results from One-Dimensional Model . . . . . . . . . . . 273.2 Analytical Model . . . . . . . . . . . . . . . . . . . . . . . . . . 303.3 Consequent Limitations . . . . . . . . . . . . . . . . . . . . . . . 333.3.1 Mechanical stability limitation . . . . . . . . . . . . . . . 333.3.2 Positive flow and U-tubing . . . . . . . . . . . . . . . . . 344 Two-Dimensional Model . . . . . . . . . . . . . . . . . . . . . . . . . 364.1 Model Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . 364.1.1 Non-Dimensionalisation . . . . . . . . . . . . . . . . . . 384.1.2 Hele-Shaw Displacement Model in Eccentric Annulus . . 444.1.3 Closure model for τw . . . . . . . . . . . . . . . . . . . . 464.2 Computational Detail . . . . . . . . . . . . . . . . . . . . . . . . 464.2.1 Slice Model Solution . . . . . . . . . . . . . . . . . . . . 465 Results from Two-Dimensional Model . . . . . . . . . . . . . . . . . 486 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 646.1 Summary of Foamed Cementing . . . . . . . . . . . . . . . . . . 646.2 Limitations of the Model & Future Directions . . . . . . . . . . . 66Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68A Numerical Simulation Parameters . . . . . . . . . . . . . . . . . . . 73viList of FiguresFigure 1.1 Canadian primary energy consumption by type in 2016 . . . . 2Figure 1.2 Schematic of the different stages of the primary cementing pro-cess . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4Figure 2.1 20% quality foamed cement . . . . . . . . . . . . . . . . . . 10Figure 2.2 A schematic of foam viscosity versus foam quality . . . . . . 13Figure 2.3 Dimensionless consistency and yield stress of foam as a func-tion of foam quality . . . . . . . . . . . . . . . . . . . . . . . 14Figure 2.4 Comparison of different 1D foam flow hydraulics . . . . . . . 19Figure 3.1 1D simulation of displacement of example (i) . . . . . . . . . 28Figure 3.2 1D simulation of displacement of example (ii) . . . . . . . . . 29Figure 3.3 1D simulation of displacement of example (iii) . . . . . . . . 29Figure 3.4 Results from analytical displacement simulation of example (i) 32Figure 3.5 Fluid densities plotted from one-dimensional simulation of ex-ample (i) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32Figure 3.6 Density of the foam at the interface . . . . . . . . . . . . . . 34Figure 4.1 Schematic of the cementing geometry before and after unwrap-ping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37Figure 4.2 Schematic of the well and annular geometry . . . . . . . . . . 38Figure 4.3 Geometry of the eccentric annulus mapped to the Hele-Shaw cell 44Figure 5.1 2D simulation of displacement of example (i) with no eccentricity 49Figure 5.2 2D simulation of displacement of example (i) with e = 0.1 . . 52viiFigure 5.3 2D simulation of displacement of example (i) with e = 0.3 . . 54Figure 5.4 2D simulation of displacement of example (i) with e = 0.5 . . 55Figure 5.5 Incompressible displacement results . . . . . . . . . . . . . . 57Figure 5.6 2D simulation of displacement of example (ii) with no eccen-tricity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58Figure 5.7 2D simulation of displacement of example (iii) with no eccen-tricity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Figure 5.8 2D simulation of displacement of example (i) with pˆa = 20 atm 62Figure 5.9 2D simulation of displacement of example (i) with pˆa = 40 atm 63viiiNomenclatureα Gas volume fractionq Foam qualityQl Volumetric flowrate of the liquidQg Volumetric flowrate of the gasρl Liquid densityρg Gas densityρc Density of the pure cement slurryM Molar mass of the gasR Universal gas constantT TemperatureYg Gas mass fractionM˙l Mass flowrate of the liquidM˙g Mass flowrate of the gasM˙ Total flowrateCa Capillary numberσ Interfacial tensionγ˙ Shear rateτY Yield stressκ Consistency factorixn Power-law indexRe Reynolds numberc j Volume fraction of fluid jYj Mass fraction of fluid jA Cross-sectional areaS Length of perimeteru Velocity vectoru Radial velocityv Azimuthal velocityw Axial velocityβ Angle of inclinationµ Fluid viscosityg Gravity of Earthτw Wall shear stressHe Hedstrom numberpc Inflow pressurepa Back pressurepbh Bottom hole pressurezc Casing surface heightza Annulus surface heightzbh Bottom hole heightr Pipe radiusri Inner radius of the annulusro Outer radius of the annulusH Channel half widthe Eccentricity of the annulusxt TimeFr Froude numberΨM Mass streamfunctionBT Buoyancy termsCT Curvature termsST Stress termsIT Inertial termsxiAcknowledgmentsFirstly, I would like to express my deepest appreciation to my advisor, Prof. IanFrigaard for his support and for his great ideas. I have been extremely lucky to havea supervisor who cared so much about my research work and helped me shape myprofessional career.I would also like to thank my wonderful lab-mates in Complex-fluids lab.Thanks for the fun and support. I greatly look forward to having all of you asmy life-time friends in the years ahead. In particular, I will forever be grateful toAmir Maleki for his great suggestions and for his patience with my endless streamof questions. Words cannot express my appreciation to you for your assistance.I am especially grateful to my family. My mom and my brother for all the love,comfort, and support they have given me throughout my time in graduate school. Icould not have done it without them.This research has been supported financially by the British Columbia Oil andGas Commission, by BC OGRIS through project EI-2016-10 and partly by NSERCand Schlumberger through CRD project 444985-12.xiiDedicationTo my beloved brother, Navid, who has been a never-ending source of love, en-couragement and support. I am truly thankful for having you in my life.xiiiChapter 1Introduction1.1 The Oil and Gas Industry in CanadaCanada with over 500,000 drilled wells is among the top ten producers of naturalgas and oil in the world. More than 400,000 people are involved in Canada’s oiland gas industry directly and indirectly. Moreover, this industry has a considerablecontribution in taxes and land sales, which directly affect government revenue andthe Canadian economy. This is especially true in Western Canada, even thoughthe large drop in oil price has caused a decline in both the market and oil & gasindustry, since late 2014.Energy consumption in Canada per capita is amongst the highest in the worldand this partly due to the location and cold weather. As shown in Figure 1.1, in2016 more than half of primary and secondary energy consumption came from oiland natural gas. The distribution of energy consumption shown in Figure 1.1 hasbeen almost constant during the last two decades, with an slight increase in naturalgas, renewable and alternate sources of energy; see [1]. As a result, efficiencies andimprovements in the upstream petroleum industry in Canada have the potential fora high impact on both total energy consumption and its environmental impact.Primary cementing of oil & gas wells plays an important role in the oil and gasindustry. Since primary cementing prevents subsurface hydrocarbons from leak-ing to surface, increases well production, and also seals the outside of the casing,it is a necessary step in all oil & gas well construction, production and eventual1Figure 1.1: Canadian primary energy consumption by type in 2016. Renew-able energy includes hydro, biomass, solar and wind energy. Graph isreplotted from [1].abandonment.About one third of the 500,000 hydrocarbon wells in Canada are not producing.Some are fully abandoned (i.e. sealed). Many of them are shut-in and holding pres-sure at the surface, or openly leaking. In terms of well leakage, gas leakage is morecommon and noticeable. It can cause health, safety and environmental problems.Failure of primary cementing of oil and gas wells is the number one reason forleakage of oil and gas wells. This suggests that there are still many aspects of ce-ment placement that need more research attention to improve. Different techniquesare used in primary cementing, due to different opinions on their effectiveness fordifferent situations [15]. Foamed cementing, the subject of this thesis, is one ofthese techniques.1.2 Primary CementingAn important stage of oil well construction is primary cementing in which cementis pumped to fill the annular space outside of a steel casing (or liner). Primary2cementing, which has been used for about 100 years, is a common method forplacement of cement and removing drilling mud. In placing the cement, ideally weremove all the drilling mud. The main purposes of primary cementing are: (i) toprevent well collapse by mechanically supporting the wellbore and (ii) to providea hydraulic seal in the annulus to separate the fluid-bearing zones and prevent anyuncontrolled flow behind casing. Failure of proper sealing reduces the productivityof the well and causes environmentally damaging leakage [30].The procedure of primary cementing is simply described as follows. Whendrilling a section of the well is finished, the well is filled with drilling mud. Thedrill pipe is removed from the well, the casing is inserted into the well and then aseries of fluids are pumped inside the casing. These fluids push the drilling muddown inside the casing and then, upon reaching the bottom of well, the fluids returnup in the annular section between the casing and the wellbore. The last fluid whichis a cement slurry, stays in the annular section, waiting to solidify; see [7, 30].In terms of dimensions, the diameter of a typical wellbore can start from 30-50cm and end at 10-20 cm, with a gap width of roughly 2 cm between outside of thecasing and the wellbore. The length of each section of casing (or liner) is about 10m and a number of casings are joined together, according to the required length ofthe well. The length of each cemented section is typically 100-500 m. In order tocentralize the casing in the well, at different locations centralizers are fitted to theoutside of the casing. However, it is still very common for the casing to form aneccentric annular space with the wellbore wall.3Drill new stage of well Remove drillpipe Insert steel casing Pump spacer fluid Pump lead & tail slurry Displace mud in annulus End of operation Figure 1.2: Schematic of the primary cementing process, showing the differ-ent stages in cementing a new casing (left to right).The rheology and density of the sequence of the pumped fluids should be en-gineered such that the drilling mud is entirely removed from the casing, leaving nochannels of mud. There are some general guidelines to improve mud displacementand cementing quality. Most notably, in laminar flows there should be hierarchyof density and viscosity. This means that each fluid must be denser and generatemore frictional pressure than its predecessors, [7]. Normally, the cement slurryand drilling fluid are prevented from direct contact during the primary cement-ing procedure since these two fluids are chemically incompatible. Incompatibilityand consequent contamination can cause the mixture of cement and mud to havea shorter thickening time than pure slurry, be more viscous, and/or to acceleratethe hydration. Alternatively, the contaminated mixture might not harden at all.Therefore the cement slurry and drilling mud should be physically separated withanother fluid (preflush or spacer) which is typically pumped ahead of the cement. Aspacer fluid can form the desired interface for optimized mud removal and superiorcement placement [22]. Other guidelines include avoiding eccentricity as muchas possible, mud conditioning, cement conditioning and movement of the casing.4Although these guidelines are based on some physical arguments, fundamental un-derstanding is still incomplete, often leading to poor cementing results [6, 29].In addition to geometrical complexities (e.g. skewness, eccentricity, tapering ),the fluids involved in cementing possess different physical and rheological proper-ties. Most importantly, some of them exhibit yield stress behavior; i.e. they do notflow unless the exerted shear stress exceeds a critical limit called the yield stress.These complicated effects need to be considered in any model intended to predictmud removal and cement placement.1.3 Why Foamed CementAlthough designing foamed cement is more complicated and needs specializedequipment, there are many perceived advantages in using foamed cement in pri-mary cementing; see [24].• Lightweight: Foam density can be controlled by gas content which governsthe porosity of the solidified cement. The main idea of considering foamedcements as an alternative for conventional cement, is to have lightweightcements and be able to have a control over the density of cement whichaffects bottomhole pressure.• Improved ductility: One of the biggest problems facing primary cementingis the fragility of non-foamed cement. Experiments indicate that bubbly ce-ment with 20-35% gas volume fraction in the mixture, can be up to oneorder of magnitude more flexible, compared to conventional cement. Usingbubbly cement, which has a better flexibility, can prevent cracking and fail-ure in primary cementing and reduces formation damages. However, thereare drawbacks: foamed cement with large gas volume fraction becomes toospongy to isolate fluid zones, while low gas fractions can be fragile. Thereis also an increased initial cost of nitrified cement: more than conventionalcement.• High strength to density ratio: Test results show that foamed cements even atvery low densities are sufficiently strong for drilling purposes in petroleumindustry.5In addition to the properties mentioned above, some researchers working onfoamed cement systems (and also some companies producing foamed cement)claim that the following properties also make foamed cement an ideal choice.• Superior zonal isolation: Kopp et al. [24] investigated zonal isolation ofconventional and foamed cement. They remarked that light weight, tensilestrength and other mechanical properties of foamed cement make it a greatoption for zonal isolation.• Better mud displacement: Foamed cement slurry provides a great advantagein drilling processes and may have better displacement features compared toincompressible cement, due to higher viscosity [2] .• Gas migration control: Benge et al. [4] claimed that the gas phase of the foamcontinues to expand when the volume of the foam is decreasing, and thisexpansion allows the foam pressure remain relatively constant. As a result,gas migration could be adequately controlled, and it can reduce migrationchannels in hardened foamed cement.It is also important to note that, although nitrified cement is more expensive thannon-foamed cement, nitrified cement costs much less than other flexible cementreplacement systems, such as vulcanized rubber [24].1.4 Thesis Objectives and OutlineThe objectives of this thesis are to develop a model of foamed cements in primarycementing of narrow eccentric annular geometries. The purpose is the model is tobe able to analyze displacement flows and flow-related features that might eventu-ally result in defects in the primary cementing job. Although targeting the annularflow, this model needs also to account for the flow downwards inside the casing.As well as deriving the model and solving it numerically, the aim is to apply themodel to realistic process situations to gain new insight into foamed cementing.In the next chapter we review studies related to foamed cementing. A one-dimensional model, including derivation, results and limitations, is given in Chap-ter 3. In Chapter 4, we derive the dimensionless governing equations for compress-ible fluid (foam) flow in the narrow annulus. Unlike previous work on modelling6primary cementing displacement flows, in this study the density is not constant andindeed density variations are a crucial part of the model and physical behavior. InChapter 5 numerical results and example solutions of the two-dimensional annulardisplacement equations are presented. The results are discussed, conclusions of thethesis and a summary of the contributions made are presented in Chapter 6.7Chapter 2Literature reviewHere we review a range of work relevant to the use of foamed cements in primarycementing operations.2.1 What Are Foamed CementsFoamed cement is created when a gas (usually nitrogen), is injected to the purecement slurry at high pressure. In order to have a light stable foam, adequateamounts of foaming agent and stabilizer should be added to the mixture as well.The created foam, contains small separate bubbles that should not coalesce andare not interconnected [12]. All foam properties strongly depend on the volumefraction of gas α . However, apart from the density, there is significant difficulty incharacterizing these properties. At least two reasons underlie this difficulty.First, much of our understanding of foam mechanics comes from study of clas-sical systems such as soap film foams, which are characterized by very high gasfractions (α > 0.95) and bubbles separated by thin films of liquid. The bulk flowrheology of these systems is characterized as a yield stress fluid, but the source ofthe foam structure is capillarity. In contrast, although foamed cement slurries maystart at surface with high volume fractions as they are formed (e.g. 75-90%), theseare then pumped at pressure downhole, where hydrostatic pressure compresses thegas phase significantly with depth. Downhole gas fractions in foamed slurries aremore commonly in the range 5-40%, with a rapid increase in α to surface values8occurring in the first few hundred meters.Second, unlike many other aqueous foams, here the liquid phase is itself non-Newtonian. Oilfield cements are fine colloidal suspensions with yield stress andshear-thinning behavior, often modeled as Herschel-Bulkley fluids. Part of thestructure of the foamed slurry, especially at lower α , is due to the cement rheol-ogy. Thus, to properly model foamed cement rheology we need to bridge betweenlimiting models for conventional foams (α . 1) and for bubbly liquids (α ∼ 0),while also taking into account the liquid phase rheology. This thesis is howevernot about the rheology of foamed cement. As our foamed cementing model is onlypreliminary, we are concerned mainly with vertical wells. In such wells the domi-nant contribution to the pressure comes from static pressure and not the frictionalpressure (typically 10% or less of the total pressure drop).Operationally, most of the technical literature uses the quality, q, to describethe foamed slurry1:q =QˆgQˆg+ Qˆl, (2.1)where Qˆg and Qˆl denote volumetric flowrates of gas and liquid phases, respectively.Lets suppose that a sequence of K fluids are pumped, with the last fluid being thefoamed cement. Assuming that the bubbles in the foam are well mixed and the twophases do not slip or coalesce, then q = α . This is the common assumption madefor stable foamed cements. With this same mixture assumption, the density of thefoamed slurry ρˆK , is given by:ρˆK = αρˆg+(1−α)ρˆc. (2.2)The cement density ρˆc is constant, whereas the gas density ρˆg is given by its equa-tion of state. For simplicity here we take:ρˆg =pˆMˆRˆTˆ, (2.3)where Mˆ is molar mass of the gas, Rˆ is the universal gas constant, and Tˆ is the1We adopt the notation of showing dimensional parameters with a “hat” (e.g. pˆ) and dimension-less parameters without (e.g. p).9temperature. Adoption of the ideal gas law for N2 is reasonable for well depths upto ≈3 km. To improve this description the Z-compressibility factor can be used.As the foam compresses, α (q) changes along the flow path: the gas flowrateQˆg at depth is not the same as at surface. It is instead the mass fraction Yg, of gas inthe foam, which is transported with the foamed mixture. In terms of Yg, the cementdensity is given by:ρˆK =1Yg/ρˆg+(1−Yg)/ρˆc , (2.4)andYg =αρˆgαρˆg+(1−α)ρˆl(=QˆgρˆgQˆgρˆg+ Qˆl ρˆl=ˆ˙Mgˆ˙Mg+ ˆ˙Ml), (2.5)with the bracketed expressions true for a homogeneous mixture. Here ˆ˙Mg and ˆ˙Mlare the mass flowrates of gas and liquid phases.Figure 2.1: 20% quality foamed cement. Picture is taken from [24].A foamed cementing job is relatively complicated, since it is critical to controlthe exact proportions of stabilizer, foaming agent, gas (nitrogen), and base slurry(mixture of dry cement and water). For all of the steps of a foam cementing job anaccurate process control system is required in order to measure, calculate and keeptrack of injected foam quality at given temperature and pressure. Thus all injectionrates should be controlled and automated. A common foamed cementing systemcontains the following equipment:10• Cement slurry pumping device• Liquid nitrogen storage tank• Nitrogen pumping device• Injection pumps for both foaming and stabilizing agents• Foam generator [9]2.2 Rheology of Foamed CementUnderstanding foam rheology can be critical and the performance of foam is ratherdifficult to predict. The foamed cement can be thought of as suspension of dis-persed bubbles in a yield stress fluid. The rheological properties of a two phasebubbly mixture depend on different factors, such as: properties of the base fluid,the diameter of the bubbles and bubble stiffness [14]. Rheological parameters mustbe quantified through developed analytical models and/or fitted from experimentalresults.The approach that has been used in this thesis is very similar to the modeldescribed by the well known Krieger-Dougherty law for a suspension of discretesolid particles in a Newtonian fluid. Here we consider that the viscosity and yieldstress of the foam are functions of quality (i.e. mimicking solid volume fraction).However, on the surface of a bubble we expect slip, whereas we have no-slip fora solid particle surface [14]. Thus, we expect some underlying differences in thismodel. Furthermore, since cement is a non-Newtonian fluid, the foamed cement isanticipated to exhibit additional non-Newtonian rheological behavior as well.It has been indicated by Rust and Manga [36] and by Llewellin et al. [25], thatthe viscosity of a foamed fluid at low steady shear rate grows with the quality of thefoam and reduces at large shear rates. They have also emphasized the contributionof the bubble stiffness (i.e. deformation) to foam viscosity. Rust and Manga [36]have assumed that at high shear rate bubbles are lengthened in the flow directionand at small shear rate bubbles are spherical. Thus a capillary number has been in-troduced into mechanical models, to measuring the stiffness of bubbles in a bubbly11fluid. This is the ratio of the viscous stress to the interfacial stress.Ca =ηˆ ˆ˙γσˆRˆ, (2.6)where Rˆ is the radius of the bubble, ηˆ is the viscosity of the base fluid, σˆ is the in-terfacial tension between the two phases (liquid and gas), and ˆ˙γ is the applied shearrate. Testing has shown that a small capillary number implies that the suspensioncontains rigid non-deformable bubbles, while large capillary number indicates softand deformable bubbles in the fluid. Hence for foams with small capillary number,it can be concluded that the bubbles are stiff compared to the base fluid viscosity.In the case of a yield stress fluid one might want to compare directly the capillaryand yield stresses. For this, the plastic capillary number is computed as:Ca =τˆY (0)σˆRˆ(2.7)Where τˆY (0) is the yield stress of the suspending fluid.As mentioned above, a few assumption have been made in order to develop anduse these rheological models for bubbly fluid, which we summarize as follows.• Foam fluid is made of discrete and spherical bubbles.• Bubbles do not coalesce and are not interconnected.• Bubbles are relatively rigid and non-deformable.• Capillary number of the foam is relatively small (plastic capillary number issmaller than 0.1).• Gas volume fraction in the bubbly fluid is limited up to 50%, so that bubblesare not deformed by geometrical restrictions [14].Generally the viscosity of a foam could be up to several order of magnitudelarger than the viscosity of suspending fluid. Experimental results have indicatedthat the viscosity of the foam, which depends on foam quality and nitrogen injec-tion rate, rises with foam quality at low to medium foam quality. However, the12Considerable surge Considerable dropGradual riseQualityµ (apparent foam viscosity)Figure 2.2: A schematic of foam viscosity versus foam quality at fixed totalinjection velocity. Picture is replotted from [18].growth at small gas volume fraction is slight and in medium volume fraction ofgas it becomes quite dramatic. Finally, for foam, with high gas volume fraction,viscosity falls sharply with increasing foam quality [16, 18].Figure 2.2 from [18] clarifies this typical description further. Three differentsurfactants (Aquet 944, Cedepal FA-406, and Stepanform-1050) were mixed withnitrogen in order to provide different kinds of foams for this experiment. Aquet944, Cedepal FA-406, and Stepanform-1050 are anionic surfactants that are typ-ically used for drilling applications. Recognizing the underlying reasons for theviscosity changes in the three zones shown in Figure 2.2 is more complex than wecan review here. It would need an analysis of the activity and interaction of bub-bles within the foam and also information about stability of the foam during shearflow. More generally, experiments on foamed cement systems show that for qualityranges between 10 to 65%, the viscosity of foam increases significantly with thefoam quality [2].There is not, to our knowledge, any comprehensive model that accounts fullyfor foamed cement rheology. It is therefore fortunate that frictional pressure gra-dients are significantly smaller than static pressure gradients and may be ignored13for many foamed cementing scenarios. Essentially this means that the static pres-sure gradients mostly determine the pressure, hence gas volume fraction. It is onlyin wellbore applications where rheology is critical to the pressure or other flowfeatures, that we need a precise constitutive law depending on volume fraction.In this thesis, for the purposes of model closure, we adopt the following con-stitutive equations based on a micro-mechanical estimation of the behavior of sus-pensions of bubbles in yield stress fluid [14], which are in good agreement withexperimental results. The yield stress and consistency of foam, as a function of thegas volume fraction, are:κˆ(α) = κˆ(0)(5+3α5−2α )n+12 (1−α) 1−2n2 (2.8)τˆY (α) = τˆY (0)√(1−α)(5+3α)5−2α (2.9)Where κˆ(0) and τˆY (0) are the consistency factor and yield stress of the purefluid (slurry), respectively. It is notable that (2.8) and (2.9) have been proposed asvalid for cases when the quality of the foam is less than 50%. We do not includeany variation in power-law index with quality.Figure 2.3: Dimensionless consistency (for n= 1) and yield stress of foam asa function of the gas volume fraction. Both graphs replotted from [14].Equation (2.8) with n = 1 and (2.9) are plotted in Figure 2.3. As observed,14yield stress is not greatly affected by the existence of bubbles in the fluid, whileconsistency factor increases significantly with increasing quality [14].2.3 Operational Methods in Foamed Cementing JobTwo modes of operational control are common in foamed cement placement.Mode 1: Constant nitrogen to base slurry flowrate ratio. Nitrogen is pumped atconstant volumetric flowrate (relative to the base slurry volumetric flowrate)on surface where a set temperature and pressure are applied. The pressureand temperature fix the gas density and hence the inflow mass fraction Yg isspecified. There are two ways of applying this.Mode 1A: the pressure is specified at the top of the casing (inflow pressure).Mode 1B: the pressure is specified at the top of the annulus (back pressure).We will use both the above in our model. Although mode 1A is easier touse mathematically, operationally a back pressure is often applied at the topof the annulus to prevent quick depressurization of the foam. Mode 1 istypically elected when cementing a conductor or surface casing for deepoffshore wells. A drawback of mode 1 is that large density gradient canresult in the annulus.Mode 2: Constant density: Density in the annulus at end of placement is keptmore or less constant (e.g. ≈ 100 kg/m3 tolerance) by adjusting the foamquality on surface. In this situation depending on the well, we can end upwith more than 20 pump stages with different nitrogen flowrates. The morenumber of stages of injecting nitrogen along the well, the higher homogene-ity of the foam density. Here back pressure is either not applied or is heldconstant. The problem is to design the pump schedule, say in terms of stepsin Yg at fixed surface temperature and pressure. This mode of operation isencountered in land wells, in the United States and Canada [28]. We havenot implemented this mode in our models.152.4 One-Dimensional Hydraulic Models of Foam FlowThe primary cementing operation using foamed cements is conceptually similar tothat using an incompressible cement slurry, in that a sequence of fluids is pumpeddown inside the casing to bottom-hole, through the float collar/shoe and returns upthe annulus. As with conventional cementing, the foamed slurry is preceded by oneor more preflushes (e.g. wash, spacer) according to the job specifics. Additionalcomplications are physico-chemical, mechanical and equipment related.(i) To ensure foam stability it is critical to control the exact proportions of stabi-lizer, foaming agent, gas (nitrogen), and base slurry (mixture of dry cementand water). The preflushes may also need to be designed for chemical com-patibility, avoiding degradation of the foam during setting as a degree ofmixing inevitably occurs.(ii) Fluid-mechanical complications arise from compressibility of the foamed slurryand from the complex dependence of the rheological properties of the slurryon the gas fraction.(iii) In terms of equipment, as discussed in §2.1, specific foam cementing systemis required in order to generate the foam and control the parameters duringthe foam cementing job. Frequently, the foamed cement job is pumped in amanaged pressure environment, with back-pressure held on the annulus, soas to control compression (i.e. density) along the flow path.The general method widely used in drilling and petroleum industry, for foamdrilling hydraulics is a pressure traverse calculation technique. In this method, thecalculation is started from the top of the casing (inlet) where the pressure is knownor assumed to be given and the entire path (pipe and annulus) is simplified as aone dimensional region. The pressure drop over the length at given total massflowrate is calculated by integrating the total pressure gradient. The total pressuregradient is the sum of frictional pressure gradient, hydrostatic pressure gradient,and acceleration/deceleration pressure gradient contributions (often neglected insteady flow calculations). The calculation continues downwards inside the pipe(casing) and then upwards along the annulus until we get to the top of the annulus(outlet). At this point, if the pressure condition is not satisfied (unequal to the16given back pressure at outflow), the whole calculation must be repeated with a newassumed inlet pressure. Since the gas phase of the foam mixture is compressible,the foam density and quality should be updated at each step with the new pressure[11].We now review a few different foam drilling hydraulics models. The aim ofeach and every drilling hydraulic model is to find a method for calculating thepressure drop along the flow path (over the length of well/casing) [31].• Beyer et al. [5] studied foam flow behavior by focusing on wall slip. Theslip component of velocity was considered to be additive to the fluid veloc-ity: u = u f luid +uslip. They assumed that the foam behaves like a Binghamplastic fluid and derived an equation for calculating the slip velocity as afunction of foam quality and wall shear stress. Using the formula of theBingham fluid velocity, they found a relationship between viscosity of thefoam and its quality. Ultimately, they presented a model for calculating fric-tional pressure drop as a function of velocity (slip velocity + fluid velocity),foam quality, and the diameter of the pipe.• Two years after Beyer’s model was published, Blauer et al. [8] started towork on predicting friction losses of a foam flowing in a pipe for differentflow regimes. This group also considered foam as a Bingham fluid. Theyoffered correlations for calculating the density and viscosity of the foam us-ing foam quality and calculated the Reynolds number and Fanning frictionfactor as a function of foam viscosity, density, flow rate and diameter of thepipe. Finally, they presented a model for calculating the frictional pressurelosses of laminar flow in a pipe, by modifying the Buckingham-Reiner equa-tion. The contribution of the gas phase to the density was neglected in thismodel.• Later in 1983, Sanghani & Ikoku [37] completed Blauer’s work by modellingfoam as a pseudoplastic fluid. He used experimental data from pipe flowof a pseudoplastic fluid to fit a model and then calculated the consistencyand power-law index of foam as a function of the gas volume fraction. Theother difference between his work with the previous model is that he did notneglect the effect of the gas phase on the foam density.17• In 1986 Reidenbach et al. [35] derived an empirical model for calculatingrheological parameters of carbon dioxide and nitrogen foams. Reidenbachet al. defined the apparent viscosity of foam as a function of pipe radiusand the rheological parameters of foam (yield stress, consistency factor, andpower-law index). Then he suggested to use the apparent viscosity in theregular pressure drop formula.• In the early 1990s Valko and Economides [41] proposed a new model forfrictional pressure gradient calculations of isothermal steady-state flow ofcompressible fluid (foam) in a pipe. They introduced a variable called thespecific volume expansion ratio, which was defined as the ratio of the liquidphase (base slurry) density to the foam density, and they used this variable toscale all density dependent parameters. Hence a new definition was offeredfor calculating volume-equalized Reynolds number, Fanning friction factor,and frictional pressure drop. Unlike other studies that had used foam qualityin their models, the key parameter in this model was the specific volumeexpansion ratio.• In late 1990s Gardiner et al. [19] continued working on the model developedby Valko and Economides. They considered foam as a psuedoplastic fluid.A modified Hagen-Poiseulle equation was derived. The effect of slip wasalso considered in the model.• Ozbayoglu et al. [31] performed a comparative study of the one-dimensionalfoam flow models.They developed their own experimental setup and also acomputer program, in order to compare the existing models with the exper-imental results. Comparison of the above models with experimental data,shows that the values of frictional pressure gradient calculated using thesemodels can be very different to actual values. Figure 2.4 shows an exam-ple comparison of the pressure drops obtained by using different foam flowmodels for a 70% gas volume fraction foamed fluid in a 10 cm diameter pipe(from [31]).18SanghaniBlauerBeyerValkoGardinerQ (m^3/s)Pressure drop (kg/m^2)0 0.005 0.01 0.015 0.02120001000080006000400020000Figure 2.4: Comparison of different one-dimensional models for foam flow[31].• A relatively recent foam drilling hydraulics and rheology model is that ofChen et al. [11]. Chen et al. studied foam pressure drop and rheology forlow quality foams flowing in horizontal, vertical or an inclined wellbores.They considered a power-law model for foam rheology. In addition, the con-sistency and power-law index were modeled as a function of foam quality,liquid phase properties (viscosity), and some empirical parameters.In conclusion, we see both an evolution in these 1D models and a wide range of pre-dictions (Figure 2.4). The differences in predictions from these one-dimensionalmodels can only come from different hydraulic closure models for the foam den-sity and rheology, (or the neglect of terms in the momentum balance, e.g. acceler-ation/deceleration).192.5 Flow in the AnnulusIn the context of annular flow, there have been studies of (single phase) non-Newtonian fluid flows since the 1960s, with the aim of allowing for hydraulicsstyle calculations. Notable here is the work of Hanks and co-workers (e.g. [23])who considered laminar, transitional and turbulent regimes in concentric annuli. Inthe cementing context, since the annulus is narrow it has been common to use “slotapproximations”, i.e. modelling the flow locally as a plane channel and neglectingcurvature. A review of hydraulics approximations can be found in [30] and morerecent approaches are reviewed in [26].When we move to flow in an eccentric annulus planar flow approximationsremain valid for as long as the annulus is narrow (as common in cementing) andare still widely used. The first reliable computations of full annular flows of yieldstress fluids were conducted in the early 1990s; see [39, 42]. The narrow gapapproximation for annular flows of Newtonian fluids was extended by de Pina etal. [13] to any ratio of inner to outer radius. Later in 2009 Gomes and Carvalho [20,21] extended dePina’s model to non-Newtonian fluids, by defining a correlation forcalculating the viscosity of the non-Newtonian fluid and using this viscosity in theequation derived to predict the Newtonian fluid pressure gradient.2.5.1 Two-Dimensional Flow of Incompressible Fluid in AnnulusIn moving to look at displacement flows in the annulus, we review only the two-dimensional models developed to date. Although a few authors have conducted 3Dsimulation studies using CFD tools, these are limited to study of a small part of theannulus and not feasible over the wellbore scale.Some of the earliest work is due to Tehrani et al. [40], who analyzed laminardisplacement flows in an annulus from two points of view: experimental and the-oretical, finding reasonable agreement. They performed experiments in a lab-scalefully inclinable annulus and simulated the same experiments. They reported resultsthat explored variations in: inclination angle, spacer design and density contrast,examining the effects on the mud removal process.One of the proposed annular displacement models that has been used in thisthesis (and modified for compressible fluid displacement) is Bittleston et al.’s model [7].20A number of assumptions have been made in this model in order to simplify thethree-dimensional annular flow to a two-dimensional flow. For instance, by assum-ing fast radial diffusion and a narrow gap between the casing and the hole, theysimplified the problem which resulted in a series of two-dimensional partial dif-ferential equations, similar to those derived later. Pelipenko and Frigaard [32, 33]performed an extensive analysis of the model derived in [7]. They derived designconditions on the rheology and density such that the displacement flow in an ec-centric annulus could still be stable. When these conditions could not be satisfied,the flow becomes unsteady, with the interface advancing rapidly along the wideside of the annulus. In strongly eccentric geometries, the narrow side fluid fails tomove at all and a narrow side mud channel can be formed.Carrasco-Teja et al. [10] studied laminar flow in horizontal eccentric annuli,again using the model of [7]. They showed that most of the time, there are steady-state traveling wave solution for a displacement flow. And this may also happento the cases when density difference between the fluids is large or even when theannulus is extremely eccentric. They discovered that when the displacing fluid isheavier than the displaced fluid, steady-state conditions are easier to be satisfied,compared to the case that light fluid displaces heavy fluid. This is because eccen-tricity of the annulus opposes buoyancy-driven slumping in the first case. Theyfound different conditions under which the displacement flow could be steady.Very recently, Maleki and Frigaard [27] have extended the same modelling ap-proach to turbulent displacement flows, exploiting the narrow gap of the annulusto reduce the momentum equations locally to a 2D turbulent shear flow. After sim-ilar calculations to [7] this again becomes an elliptic PDE for the streamfunction.They present a number of interesting examples with turbulent and mixed regimes.A main difference with the laminar flows is that dispersion/diffusion between thefluids is now a significant effect and must be modelled.The above cited works have been developed largely at UBC and collaborativelywith Schlumberger. There are a large number of papers in the oilfield literature thatpresent examples of the usage of these models for field case studies. In recent yearsa number of other groups have developed broadly similar models. The work ofCarvalho and co-workers, plus collaborators in Petrobras is documented in [20, 21]and has some nice experimental comparisons in [3].21Other proprietary models exist within other companies but the models are notexplained in the open literature, nor openly validated. Some models are mar-keted as “3D”, but in fact only reconstruct the velocity variation across the gap,from a simpler flow. This can of course be done with any of these 2D gap-averaged models, including that developed later in this thesis. Others models,e.g. [20, 21], claim additional validity for wider gap geometries. This is also onlypartly true. Most primary cementing flows are inertial with significant Reynoldsnumbers (e.g. Re> 50). Without the narrowness of the annulus, inertial effects willlikely be present in these flows limiting the benefits of improved friction factor clo-sures.22Chapter 3One-Dimensional Version of aGeneral Model3.1 One-Dimensional ModelWe now develop a one-dimensional hydraulic model to describe transport of thefluids along the flowpath. We assume that a sequence of K fluids is pumped intothe casing, each characterized as a Herschel-Bulkley fluid. For simplicity only thelast stage (fluid K) is a foamed slurry. The flow path is parameterized by zˆ, fromzˆ = zˆc = 0 at the inflow at surface, down the casing to zˆbh at bottom hole, and re-turning upwards along the annulus to zˆa at the exit of the annulus. To simplify, wewill only consider onshore wells and assume that zˆc and zˆa are at the same heighton the surface; hence zˆa− zˆc = 2zˆbh. To simplify the model we neglect temperaturevariation along the flow path. Wellbore heat transfer is usually convectively domi-nated, so this approximation degenerates for deep wells where conduction has timeto act. Thus, the assumption is that the nitrogen is injected at a specified surfacepressure pˆc (and temperature), at the top of the casing.The fluids are modelled as a locally homogeneous mixture, with c j representingthe volume fraction of fluid j. The mixture moves with a single velocity uˆ, (no slip23between the constituents or phases). Mass conservation is modeled as:∂∂ tˆ(ρˆ jc j)+ ∇ˆ · (ρˆ jc juˆ) = 0, j = 1,2, ...,K, (3.1)Where ρˆ j and c j are density and volume fraction of fluid j, respectively. Theseequations may also be written in terms of mass fraction Yj of fluid j, as follows.∂∂ tˆ(ρˆYj)+ ∇ˆ · (ρˆYjuˆ) = 0, j = 1,2, ...,K. (3.2)whereρˆ =K∑j=1ρˆ jc j =1K∑j=1Yjρˆ j. (3.3)On summing (3.2) over j we find:∂ ρˆ∂ tˆ+ ∇ˆ · (ρˆuˆ) = 0, (3.4)representing mass conservation of the total mixture. On subtracting (3.4) from(3.2) we find that∂∂ tˆYj + uˆ · ∇ˆYj = 0, (3.5)i.e. each mass fraction is advected with the flow. This may be repeated specificallyfor the slurry mass conservation equation, with density ρˆK , to show that:∂∂ tˆYg+ uˆ · ∇ˆYg = 0, (3.6)for the gas mass fraction.The cross-sectional area, at distance zˆ along the flow path, is denoted Aˆ(zˆ).In a one-dimensional model we assume that the mass (or volume) fractions varyprimarily in the zˆ-direction. On averaging across the flow path, we find:∂∂ tˆ(AˆρˆYj)+∂∂ zˆ(AˆρˆYj ˆ¯w) = 0, j = 1,2, ...,K, (3.7)24or advectively:∂∂ tˆYj + ˆ¯w∂∂ zˆYj = 0, j = 1,2, ...,K, and g (3.8)and on summing:∂∂ tˆ(Aˆρˆ)+∂∂ zˆ(Aˆρˆ ˆ¯w) = 0. (3.9)Here ˆ¯w is the averaged velocity in the zˆ-direction. We assume that the main changesin mixture density are due to advection of the mixture and approximate (3.9) by∂∂ zˆ(Aˆρˆ ˆ¯w) = 0. (3.10)Thus, the total mass flow rate across each cross-section is constant in our hydraulicmodel, given by the mass flow rate at surface: ˆ˙M = ˆ˙Mg+ ˆ˙Ml .For given pressure and mass fractions, the density is determined and hence themean axial velocity is calculated from the mass flow rate:ˆ¯w =ˆ˙MρˆAˆ. (3.11)The pressure is calculated from a steady one-dimensional averaged axial momen-tum equation:− ∂∂ zˆ(Aˆpˆ)− Sˆτˆw+ Aˆρˆ gˆcosβ = 0, (3.12)where β is the angle of inclination, between the zˆ-axis and vertical direction. Sˆ(zˆ)is the length of (wetted) perimeter at zˆ and τˆw is the wall shear stress. Equation(3.12) assumes a steady flow with slow axial variation in both the mean velocityand the axial stresses. The wall shear stress is expressed via the closures explainedin §3.1.1: primarily it is a function of the local flow velocity, mixture density and(within the foamed slurry), also the gas fraction and hence pressure.Although the frictional pressures (wall shear stress terms) introduce the addi-tional complexity of coupling the momentum balance with pressure, velocity andgas fraction, it is worth noting that the frictional contribution to the total pres-sure drop is generally small, (. 10% of the static pressure). Thus, as an initial25approximation, for many wells these terms may be neglected or evaluated in asemi-implicit way when computing the pressure. However, these terms may be-come important in horizontal wells and in all wells in situations of rapid expansionand/or unsteady flow (e.g. U-tubing, between casing and annulus if no check-valveis fitted), so our modelling approach is preliminary.3.1.1 Wall Shear StressBelow we shall want to evaluate the wall shear stress in the pipe and annulus, as afunction of the rheological parameters. Following Maleki & Frigaard’s model [26]for laminar pipe and plane channel flows, the Reynolds number for Herschel-Bulkley fluids is:Rep =8ρˆWˆ02κp( 4ˆ¯wrˆ )n, pipe flow,6ρˆWˆ02κp( 6ˆ¯w2Hˆ)n, channel flow.(3.13)Where Wˆ0 is the mean velocity in a pipe of radius rˆ or in a channel of width 2Hˆ,and κˆp is defined as follows:κˆp =κˆ(3n+14n )n, pipe flow,κˆ(2n+13n )n, channel flow.(3.14)The Hedstrom number and wall shear stress are:He =τˆY ( ρˆn(2rˆ)2nκˆ2p)12−n , pipe flow,τˆY ( ρˆn(2Hˆ)2nκˆ2p)12−n , channel flow.(3.15)τˆw =Hw( ρˆn(2rˆ)2nκˆ2p)12−n, pipe flow,Hw( ρˆn(2Hˆ)2nκˆ2p)12−n, channel flow.(3.16)26The procedure of calculating Hw from Rep is fully described in [26].3.1.2 Results from One-Dimensional ModelAlthough the two modes described in §2.3 are different operationally, the compu-tational method of solving the 1D model is similar. Mode 2 is however complex todesign the pump schedule and consequently we focus on mode 1.At any time step, we assume that we specify at zˆc: an inflow pressure Pˆc(tˆ),inflow mass fractions of each constituent, Yc, j(tˆ), and the total mass flow rate ˆ˙M(tˆ).These are calculated from the specified flow rate ratios and other pump schedule.Assuming the Yj are known at timestep n, the momentum balance is integratedin the positive zˆ-direction, using a modified Euler finite difference scheme. Thiscomputes the pressure, density and mean velocity (from (3.11)), everywhere alongthe flow path at timestep n. In the case that the pressure is fixed at the outflow(back pressure imposed), we iteratively change the inflow pressure Pˆc(tˆ) to satisfythe outflow condition. The transport equations (3.8) are hyperbolic and are solvedusing an FCT scheme [43], to give the new Yj at timestep n+1.We now present 3 illustrative examples, which we will follow throughout thiswork. For each of these we consider only mode 1 operational/boundary condi-tions. The well geometry is fixed as a 500m deep vertical (β = 0) borehole, withcasing inner, casing outer and hole diameters: 18.42, 20.32, 25.40 cm, respectively(i.e. 7.25, 8 & 10 inches). The pipe is initially full of a drilling mud (displacedfluid) with properties: ρˆ1 = 1400 kg/m3, κˆ1 = 0.02 Pa.sn1 , n1 = 1, τˆY,1 = 4 Pa. Themud is displaced by a Nitrogen-foamed cement (displacing fluid). The propertiesof the pure cement slurry are: ρˆ2 = 1800 kg/m3, κˆ2 = 0.04 Pa.sn2 , n2 = 1, τˆY,2 = 8Pa. The inlet temperature is fixed at 300K. and the following conditions are usedat surface.Example (i): pˆc = 20 atm, Qˆg = 120 l/min, Qˆl = 180 l/min.Example (ii): pˆc = 40 atm, Qˆg = 120 l/min, Qˆl = 180 l/minExample (iii): pˆc = 20 atm, Qˆg = 180 l/min, Qˆl = 120 l/minThe focus of these examples is on the effects of density and flow rate during dis-placement. The mud rheology is not extreme and provided the annulus were ad-27equately centralized we would not expect any difficulty to displace this well. Al-though the fluid properties are realistic, we have restricted to 2 fluids (mud andcement), which is not realistic.The results are shown in Figs. 3.1 - 3.3. Each figure plots density, quality,velocity and pressure in the form of spatiotemporal colourmaps showing evolutionas the cement-mud interface travels down the casing and back up the annulus. Weobserve the faster speeds of the interface in the annulus, due primarily to smallercross-sectional area (see the mud). The foam velocity decreases in the casing withdepth and increases in the annulus, as it decompresses. We see a slight asymmetryin the pressure with respect to distance from bottom-hole, at both start and end ofthe examples (when the flow loop is full of a single fluid). This asymmetry is dueto frictional pressure losses, which are negligible here.Figure 3.1: One-dimensional simulation of displacement of Example (i).28Figure 3.2: One-dimensional simulation of displacement of Example (ii).Figure 3.3: One-dimensional simulation of displacement of Example (iii).29Comparing examples (i) and (ii), both have the same surface quality at zˆc, butthe entire system is operating at increased pressure in example (ii). Thus, the ef-fects of pressure variations are generally less pronounced in example (ii): frictionalpressure effects are smaller and the variation in quality, foam density and velocityare smaller.Comparing examples (i) and (iii), we have increased the inflow quality from40% to 60%. With the relatively short well the density in the foam remains belowthat of the mud throughout the simulation, the interface being density unstablethroughout the annulus. Frictional pressure effects are comparable to example (i),except close to the top of the annulus. Here the higher quality (and the smallerarea) leads to high velocities.3.2 Analytical ModelIf the frictional pressures are neglected, we arrive at a model that can be integratedanalytically. This is advantageous for both analysis and for identifying limitationson the foamed cementing process; see e.g. [28]. Since we neglect the frictionalterms, this hydrostatic model is equally valid for either inside the casing or in theannulus. To simplify further, we also assume a uniform cross-sectional area in eachgeometry (straightforwardly extended to piecewise constant wellbore geometries)at constant inclination β . We give only the solution for K = 2 fluid stages, (buteasily extended to K > 2).First, lets assume that the foamed slurry (fluid 2) occupies [zˆc, zˆi) with thedrilling mud occupying [zˆi, zˆa], and that the pressure pˆc is imposed. For any zˆ ∈[zˆc, zˆi), we find the pressure by integrating the momentum equations, rearranged asfollows:d pˆdzˆ= ρ2gˆcosβ , ⇒ dpˆρ2 = gˆcosβdzˆ. (3.17)The gas mass fraction is fixed in the foamed slurry, and using (2.4) we find:YgRˆTˆMˆln(pˆ(zˆ)/ pˆc)+1−Ygρˆc(pˆ(zˆ)− pˆc) ={gˆcosβ (zˆ− zˆc), zˆ≤ zˆbh,gˆcosβ (2zˆbh− zˆ− zˆc), zˆ > zˆbh.(3.18)This implicit relation is inverted numerically at each zˆ, when needed, or simply at30the interface position zˆi. The right-hand side of (3.18) changes at zˆbh due to thereversal in direction of gravitational acceleration as we return up the annulus. Forzˆ ∈ [zˆi, zˆa], the density is constant andpˆ(zˆ) = pˆi+ρˆ1gˆcosβ (zˆ− zˆi), zˆ≤ zˆbh,ρˆ1gˆcosβ (2zˆbh− zˆ− zˆi), zˆ≥ zˆbh ≥ zˆiρˆ1gˆcosβ (zˆi− zˆ), zˆ≥ zˆi > zˆbh.(3.19)This gives the annulus pressure at the exit as:pˆa = pˆi+ ρˆ1gˆcosβ (zˆi− zˆa). (3.20)If zˆi = zˆa, then evidently pˆa = pˆi, and from (3.18), if 2zˆbh = zˆa− zˆc, we find: pˆi = pˆc,with the physical meaning that the foam decompresses exactly the same in theannulus as it compresses inside the casing.Figure 3.4 shows a comparison of the of numerical and analytical models forexample (i), computed previously; see Figure 3.1. At the start and finish of thesimulation, we can see the symmetry of pressure, density and quality, with respectto distance, i.e. as expected with no frictional pressures present. This is used tovalidate the numerical method used for computing the 1D model, (i.e. by settingthe frictional pressure terms to zero in the 1D model). Evidently, by comparingFigs. 3.1 & 3.4 we can see that neglect of frictional pressure is a reasonable firstapproximation for this type of well and cementing parameters.31Figure 3.4: Results from analytical displacement simulation of Example (i).Figure 3.5 shows the fluid densities plotted against zˆ and zˆi for examples (i).(a) With frictional pressure terms (b) Without frictional pressure termsFigure 3.5: Results from one-dimensional simulation of Example (i).323.3 Consequent LimitationsThere are two important conditions should be satisfied in foamed cementing job:3.3.1 Mechanical stability limitationIn a typical cementing operation the cement slurry is separated from the other flu-ids by rubberized plugs, as it is pumped down inside the casing. These minimizethe risk of mixing and contamination on the way down the well. As the pressureincreases with depth we have seen that the foam density also increases and this is amechanically stable situation. On the other hand, in the annulus the fluid stages arenot separated. As the cement-mud (or cement-spacer) interface advances along theannulus, there is also the possibility that the displacing foamed cement slurry mayhave a density below that of the fluid above it. This situation is mechanically unsta-ble and may promote buoyancy-induced mixing. Equally, depending on the foamrheology and that of the drilling mud, there may also result a negative viscosityhierarchy - also conducive to instability.We shall explore some of these possibilities in the annular displacement modeldeveloped in Chapter 4. However, for now we can give a conservative predictionof this mechanically unstable situation by using the analytical model. Note that ne-glecting the frictional pressures means that the actual pressures along the flowpathare higher (for a fixed outflow pressure) than without friction. Hence the foamedcement slurry is less dense in the analytical model than in the frictional 1D model:density unstable configurations are necessary in the analytical model for the sameto exist in the frictional model.For the same well as in examples (i)-(iii) we plot the foam density evaluatedat the interface position, for different pˆc = pˆa. The dark contour illustrates themud density. As we have seen, for the lower surface pressures imposed the foamcompresses more down hole and higher densities are achieved downhole. Thereare clearly competing objectives here: on the one hand one motivation for usingthe foamed cement is having lower densities downhole, but on the other hand thereis a risk of density driven instability in the annulus that results.33(a) Inflow quality is 0.4 (b) Inflow quality is 0.6Figure 3.6: Density of the foam at the interface, for different pˆc = pˆa.3.3.2 Positive flow and U-tubingOther constraints come from the hydraulics. Neglecting friction, the pressure at thetop of the casing pˆc must overcome the resisting pressures, which are pˆa and thenet static pressure, i.e. in order to flow in the positive direction we require:pˆc > pˆa−zˆa∫zˆcρˆ gˆcosβdzˆ (3.21)If (3.21) is not satisfied there is a net flow driven from annulus back towards thecasing. In practice a check valve may prevent this flow from occurring. If instead(3.21) is satisfied the difference in (driving-resisting) pressures is positive and willordinarily be balanced by the frictional pressures, i.e. the flow moves at a speedsuch that the frictional pressure balances.It can also happen that the static pressure terms in (3.21) are in such imbalancethat the equilibrium flow rate exceeds the capacity of the pump. In this case thefluids are in free-fall (perhaps only temporarily) and pˆc will drop. This occurs inother pumping operations where there can result a pressure imbalance (U-tubing).These situations are a limit of our simple 1D steady hydraulics model. Althoughthe model can be modified readily, retaining the acceleration terms, an additional34complication results here in that a loss of pressure decompresses the foam which isan evident risk. For the examples presented here we simply avoid regimes whereU-tubing occurs.35Chapter 4Two-Dimensional Model4.1 Model DerivationWe now turn to the main focus of this work: to develop a two-dimensional modelfor the placement of foamed cement in the annulus. Here we stick to laminar flowregimes and the model is thus an extended version of that initially proposed byBittleston et al. [7]. Much of our model derivation and simplification to a Hele-Shaw model follows the scaling arguments in [7, 27], and we are consequentlybrief in this regard. For simplicity we model only half of the annulus, assumingsymmetry at wide (φ = 0) and narrow (φ = 1) sides; see Figure 4.1. However, it isnot a strict assumption required for the modeling. In practice the same model willwork if the wide side of the annulus is at any specified φw(ξ ) (e.g. if the symmetrycondition is replaced by periodicity in φ ).A cylindrical coordinate system (rˆ,θ , ξˆ ) is used to describe the well geometry:ξˆ measures distance along the central axis of the casing rˆ = 0 which is assumed tobe inclined to the vertical with angle β (ξˆ ). In the local coordinate system the fol-lowing Navier Stokes system is satisfied by the velocity uˆ = (uˆ, vˆ, wˆ) and pressure36previous casingdrilling mudspacerLead slurryFoamed cementslurrynewly cemented casinghydraulically isolated zones(a) (b)After primary cementingFlow direction Flow direction(c) (d) Figure 4.1: Schematic of the cementing geometry before and after unwrap-ping.pˆ. The dimensional momentum equations are1:[∂∂ tˆ+ uˆ · ∇ˆ](ρˆ uˆ) = 1rˆ∂∂ rˆ[rˆτˆrˆrˆ]+1rˆ∂∂θτˆrˆθ +∂∂ ξˆτˆrˆξˆ −τˆθθrˆ− ∂ pˆ∂ rˆ+ ρˆ gˆrˆ,(4.1)[∂∂ tˆ+ uˆ · ∇ˆ](ρˆ vˆ) = 1rˆ2∂∂ rˆ[rˆ2τˆθ rˆ]+1rˆ∂∂θτˆθθ +∂∂ ξˆτˆθξˆ −1rˆ∂ pˆ∂θ+ ρˆ gˆθ (4.2)[∂∂ tˆ+ uˆ · ∇ˆ](ρˆwˆ) = 1rˆ∂∂ rˆ[rˆτˆξˆ rˆ]+1rˆ∂∂θτˆξˆ θ +∂∂ ξˆτˆξˆ ξˆ −∂ pˆ∂ ξˆ+ ρˆ gˆξˆ (4.3)1As earlier, we adopt the notation of showing dimensional parameters with a “hat” (e.g. pˆ) anddimensionless parameters without (e.g. p).37Where the components of the gravity acceleration vector are:gˆrˆ =−gˆsinβ (ξˆ )cosθ gˆθ = gˆsinβ (ξˆ )sinθ gˆξˆ =−gˆcosβ (ξˆ ) (4.4)and gˆ = 9.81 m/s2. Note, we have retained the earlier assumptions regarding mod-elling the foam as a compressible fluid with no slip between phases. Each fluidin the sequence of K fluids is characterized as a Herschel-Bulkley fluid, with thefoam rheology described as in §2.2. The fluids combine locally to form a misci-ble mixture defined through the mass fractions Yj of each constituent and thus wehave a single velocity and pressure field. Mass conservation for the total mixtureand individual components are described by (3.4)-(3.6). Temperature variationsare ignored.4.1.1 Non-DimensionalisationDimensionless Geometry Parametersˆ  ˆ ˆeˆ ˆoˆr ˆiˆrFigure 4.2: Schematic of the well and annular geometry. eˆ(ξˆ ) is the displace-ment of the two centers of the two cylinders.The local cross-section of the well is an eccentric annulus with inner radius rˆi(ξ )(outer radius of the inserted casing) and outer radius rˆo(ξ ) (equal to the radius of38the wellbore or previous casing); see Figure 4.2. We define the (local) mean radiusand mean half-gap width:rˆa(ξˆ )≡ 12 [rˆo(ξˆ )+ rˆi(ξˆ )], dˆ(ξˆ )≡12[rˆo(ξˆ )− rˆi(ξˆ )], (4.5)and then average along the flow path:rˆ∗a =1Zˆ∫ ξˆtzξˆbhrˆa(ξˆ )dξˆ , Zˆ = ξˆtz− ξˆbh (4.6)Here ξˆbh denotes the bottomhole and ξˆtz is the top of the zone of interest. Theξˆ coordinate is curvilinear along the flowpath (annulus axis), measured upwardsfrom bottom hole. Local and mean aspect ratios are defined by:δ (ξˆ ) =dˆ(ξˆ )rˆa(ξˆ ), δ ∗ =1Zˆ∫ ξˆtzξˆbhδ (ξˆ )dξˆ . (4.7)Using these parameters we may define the dimensionless (local) average radiusand eccentricity of the annulus as:ra(ξ ) =rˆa(ξˆ )rˆ∗ae(ξ ) =eˆ(ξˆ )2dˆ(ξˆ )(4.8)Dimensionless radial, azimuthal, and axial coordinates are defined by:y =rˆ− ra(ξ )rˆ∗arˆ∗aδ ∗, φ =θpi, ξ =ξˆ − ξˆbhpi rˆ∗a. (4.9)In terms of these coordinates, the external (wellbore) and internal (casing) wallsare defined by y =±H(φ ,ξ ), where:H(φ ,ξ ) =δ (ξ )ra(ξ )(1+ e(ξ )cospiφ)δ ∗, (4.10)to leading order in δ ∗/pi .The parameter δ ∗/pi is representative of the ratio of annular gap width to cir-cumference of the annulus. In a typical well the mean annular gap remains constant39while the diameter of casings decreases with depth, so δ ∗/pi may range from 0.02to 0.06 along a typical well.The narrow gap assumption, δ ∗/pi 1, is made below and used extensively tosimplify and retain only the leading order effects. The following assumptions areeither implicit in this or stated explicitly: (i) Axial variations in the cross-sectiongeometry and inclination are very slow. (ii) it is assumed that 2dˆ(ξˆ )> eˆ(ξˆ ), whichmeans the casing does not touch the well wall, i.e. the eccentricity e(ξ ) ∈ [0,1).(iii) we have assumed that the wide side of the annulus is on the upper side of thewell.Dimensionless Velocity ParametersTo develop scales we consider that a pump schedule is imposed at the top of thecasing. Via a one-dimensional model (as described in the previous section) thisdetermines the inflowing fluid(s) at bottomhole. Given that foamed slurry is com-pressible we specify the pump schedule as a mass flow rate into the well: ˆ˙Mpump(tˆ).We define a reference cross-sectional area of the annulusAˆ0 = 4piδ ∗(rˆ∗a)2, (4.11)and take as reference density the in-situ mud density ρˆ1.The dimensionless pump schedule isM˙(t) =ˆ˙Mpump(tˆ)ˆ˙M0: ˆ˙M0 = maxtˆˆ˙Mpump(tˆ), (4.12)and the velocity scale is:ˆ¯W =ˆ˙M0ρˆ1Aˆ0=maxtˆ ˆ˙Mpump(tˆ)4piδ ∗[rˆ∗a]2ρˆ1. (4.13)All azimuthal and axial velocities are scaled with ˆ¯W , and radial velocities are scaledwith ˆ¯Wδ ∗/pi . The scaled time is the axial length-scale divided by the referencevelocity:tˆ = tˆ0t tˆ0 =pi rˆ∗aˆ¯W(4.14)40Dimensionless Fluid PropertiesThe dimensionless density is simply:ρ =1K∑j=1Yjρ j: ρ j =ρˆ jρˆ1, (4.15)For j = K the foamed slurry isρK = αρg+(1−α)ρc = 1Ygρg+1−Ygρc, (4.16)where ρc = ρˆc/ρˆ1 and ρg = ρˆg/ρˆ1 is pressure dependent, following (2.3).Each individual fluid j is characterized by (τˆY, j, κˆ j,n j), indicating yield stress,consistency and power-law index of fluid j in the pumped sequence. The narrowannulus resembles locally a plane channel and so, following [27] we define a strainrate relevant to laminar flow through a plane channel:ˆ˙γ0 =ˆ¯Wδ ∗rˆ∗a, (4.17)which we use to define τˆ0:τˆ0 = maxj=1,...,K[τˆ j,Y + κˆ j ˆ˙γn j0 ]. (4.18)The stress scale is then used to scale the rheological parameters:τ j,Y =τˆ j,Yτˆ0, κk =κˆ j ˆ˙γn j0τˆ0. (4.19)Dimensionless Momentum EquationsGiven that δ0/pi  1, we may assume that the dominant components of meanvelocity will be in the (φ ,ξ )-directions, scaling approximately with ˆ¯W . Thesevelocities will give rise to viscous stresses, due to shear in the narrow gap. Using41the rheology, gap size and ˆ¯W we can develop an estimate for the size of shearstresses, say τˆ0, which we use in the scaling below. Having determined the shearstress scale τˆ0 (see §4.1.1) the relative sizes of other components of the deviatoricstress follow from the velocity and length-scales. For the pressure, we subtract offthe bottom hole pressure and a static component due to the 1st fluid, and then scalethe remaining pressure to balance the viscous stresses at leading order:pˆ = pˆbh(tˆ)+ ρˆ1g ·x+(piτˆ0/δ0)p. (4.20)We substitute into the momentum equations (4.1), (4.2) & (4.3) and dividethrough by the largest dimensional scales, to give the following.O(δ 30pi3ρˆ ˆ¯W 2τˆ0)︸ ︷︷ ︸IT= −∂ p∂y+O(δ0piρ−1Fr2)︸ ︷︷ ︸BT+O(δ 20pi2)︸ ︷︷ ︸ST+O(δ0pi)︸ ︷︷ ︸CT, (4.21)O(δ0piρˆ ˆ¯W 2τˆ0)︸ ︷︷ ︸IT= − 1ra∂ p∂φ+(ρ−1)sinβ sinpiφFr2+∂∂yτφy+O(δ 20pi2)︸ ︷︷ ︸ST+O(δ0pi)︸ ︷︷ ︸CTO(δ0piρˆ ˆ¯W 2τˆ0)︸ ︷︷ ︸IT= −∂ p∂ξ− (ρ−1)cosβFr2+∂∂yτξy+O(δ 20pi2)︸ ︷︷ ︸ST+O(δ0pi)︸ ︷︷ ︸CT(4.22)where ρ is the scaled density and Fr is the Froude number:Fr =√τˆ0ρˆ1gˆδ0rˆa,0.Our leading order model is the narrow gap limit, δ0/pi → 0, with other parametersfixed. The lower order terms that vanish in this limit are identified by the under-braces as follows: IT - inertial terms; BT - buoyancy terms; ST - stress terms; CT -terms associated with neglect of curvature. It is interesting to note that the leadingorder momentum balance is not affected by the compressibility, except through thedefinition of the density and rheology in terms of the gas fraction (quality).42Narrow gap approximationThe dimensionless continuity equation for our system is:∂ρ∂ t+∂∂y(ρu)+1ra∂∂φ(ρv)+∂∂ξ(ρw)+O(δ0pi)︸ ︷︷ ︸CT= 0 (4.23)As before, we assume that ρ changes principally via advection, ignore the timederivative above and also the curvature terms (marked CT ), taking δ0/pi → 0. Wethen integrate across the gap between y =∓H:∫ H−H1ra∂∂φ(ρv) dy+∫ H−H∂∂ξ(ρw) dy =1ra∂∂φ∫ H−Hρv dy+∂∂ξ∫ H−Hρw dy = 0.(4.24)We have used no slip to eliminate the boundary terms above. We have seen thatto leading order there is no pressure variation across the annular gap (4.21) andnow assume that the constituent mass fractions are approximated to leading orderby their gap-averaged values. Therefore, the density is simply evaluated with the(gap-averaged) pressure and mass fractions.The mass conservation equation is satisfied by using a stream function. Herehowever, we differ from [7, 27] in that our stream function ΨM is based on themass flux (not volume flux as is usual):− ∂ΨM∂ξ=∫ H−Hρv dy = 2Hρ v¯,1ra∂ΨM∂φ=∫ H−Hρw dy = 2Hρw¯ (4.25)where (v¯, w¯) are the gap-averaged azimuthal and axial velocity components. Wehave(−v¯, w¯) = 12Hρ∇aΨM (4.26)Where ∇a is the azimuthal gradient. This and the azimuthal divergence operatorare defined as:∇aX = (1ra(ξ )∂X∂φ,∂X∂ξ) ∇a ·X = 1ra(ξ )∂X∂φ+∂X∂ξScaling and averaging of the individual component mass conservation equa-43tions leads to:∂∂ tYj +(v¯, w¯).∇aYj = 0, j = 1,2, ...,K,g, (4.27)again to leading order in δ ∗/pi .Finally, taking δ ∗/pi → 0, the momentum equations become:0 = −∂ p∂y, (4.28)0 = − 1ra∂ p∂φ+(ρ−1)sinβ sinpiφFr2+∂∂yτφy (4.29)0 = −∂ p∂ξ− (ρ−1)cosβFr2+∂∂yτξy (4.30)4.1.2 Hele-Shaw Displacement Model in Eccentric AnnulusFigure 4.3: Geometry of the eccentric annulus mapped to the Hele-Shaw cell.The equations (4.28)-(4.30) can now be treated as in [7, 27], to derive a com-pressible Hele-Shaw model, which is geometrically equivalent to unwrapping thenarrow annulus into a two-dimensional channel of varying width; see Figure 4.3.To remove the y-dependency we use the condition of symmetry: (τφy,τξy) =44(0,0) at y = 0, and integrate (4.29) & (4.30):(τφy,τξy) = y(1ra∂ p∂φ− (ρ¯−1)sinβ sinpiφFr2,∂ p∂ξ+(ρ−1)cosβFr2). (4.31)As in [7, 27] we may deduce that the direction of flow is down the modified pres-sure gradient. The flow direction, say es is also given by the gap-averaged velocity,and hence stream function:es =1|∇aΨM|(−∂ΨM∂ξ,1ra∂ΨM∂φ) (4.32)By changing coordinates so that locally (s,n) are in the direction of flow, (4.31)can be reduced to a single component in the direction of es:τsy =− yH τw. (4.33)Here τw is the dimensionless wall shear stress, defined as follows:τw = H|( 1ra∂ p∂φ− (ρ¯−1)sinβ sinpiφFr2,∂ p∂ξ+(ρ¯−1)cosβFr2)| (4.34)Combining (4.32) and (4.34), we have:(−∂ΨM∂ξ,1ra∂ΨM∂φ)|∇aΨM| =−Hτw(1ra∂ p∂φ− (ρ−1)sinβ sinpiφFr2,∂ p∂ξ+(ρ−1)cosβFr2).(4.35)After cross-differentiating to eliminate the pressure, we have:∇a · [S+b] = 0, (4.36)whereS =raτw(|∇aΨM|)H|∇aΨM| ∇aΨM, b =ra(ρ¯−1)Fr2(cosβ ,sinpiφ sinβ ). (4.37)Remarkably, (4.37) is identical with that for incompressible dispacements, as in[27], except that the stream function used is the mass stream function ΨM.454.1.3 Closure model for τwHere we outline the closure relationship that defines τw(|∇aΨM|;φ ,ξ , t), locally inthe annulus, and henceS =raτw(|∇aΨM|)H|∇aΨM| ∇aΨM.we adopt the hydraulic framework introduced in [26], where the flow of a Herschel-Bulkley fluid along a narrow channel is studied in laminar, transitional and turbu-lent regimes. We are only concerned with the laminar regime. Given the shearstress scale τˆ0, the dimensional wall shear stress is τˆw = τˆ0τw. The dimensionalmean speed (averaged across the local gap) and the dimensional local annular gapare given as:Wˆ0 = ˆ¯W12ρH|∇aΨM|, 2Hˆ = 2Hδ0rˆa,0 (4.38)recall that ˆ¯W is the velocity scale for the entire annulus. We also assume that, ac-cording to the mass fractions (Yj : j = 1,2, · · ·K−1,g) at (φ ,ξ , t) we may constructthe dimensional local mixture density ρˆ and the rheological parameters: τˆY , κˆ , n.The desired closure can be obtained using the relation between Wˆ0 and τˆw whichwas given earlier in equation (3.16).4.2 Computational Detail4.2.1 Slice Model SolutionAt every time step, we solve (4.36) to update the mass stream function and subse-quently the velocity field. We then move to the next time step by updating the massfractions through solving (4.27). Note that (4.36) is not directly time-dependent:time dependency only enters indirectly through Yj or possible changes in massflow rate M˙. In the case of incompressible fluids, equation (4.36) is a nonlinearelliptic partial differential equation, which is properly formulated and solved in avariational setting due to non-differentiability associated with resolving regions ofzero flow in the annulus. For example, in [27] we have solved this using an aug-mented Lagrangian method. Here we have the additional complexity of pressure-dependency of the rheology and density. We have not so far analyzed the effect of46compressibility on the existence of solutions to (4.36), nor how this might affectmethods such as the augmented Lagrangian method. However, other pressure-dependent systems have been analyzed, e.g. [17], and we anticipate that since thefrictional effects are anyway small for most foamed cementing flows, developmentof a theoretical and computational analysis should be possible for practical flowregimes.Here however, we adopt a simpler but approximate method to solve (4.36).Due to the length of the annulus, we assume that the ξ -derivatives in (4.36) aregenerally small and neglect them, which leads to:∂∂φ[Sφ +bφ ] = 0.or in other words:Sφ +bφ = G(ξ ). (4.39)Physically, G is the modified pressure gradient and is independent of φ . Thus,on each annular slice, ξ = ξi, we iteratively change G until the net imposed massflowrate through the exit section is attained. Although this appears a convolutedprocedure, it is straightforward numerically and is robust. The advection equationfor the mass fraction (4.27) is solved using the Flux-Corrected Transport (FCT)method, [43] to capture any shock behavior. This method has been used in [27]and previously for our incompressible model simulations.47Chapter 5Results from Two-DimensionalModelTo explore different features of our model, we revisit Examples (i)-(iii), introducedearlier in §3.1.2, and compute the two-dimensional flow. Figure 5.1 shows theresults of simulation from Example (i), assuming a concentric annulus. Figure5.1b shows snapshots of the mud mass fraction Y1, at regular time intervals as thedisplacement progresses along the annulus. Superimposed on the mass fractioncontours are the streamlines of the mass streamfunction ΨM. As explained onlyhalf the annulus is considered, from wide (W) to narrow (N) side on the horizontalaxis. For the majority of the displacement the streamlines are parallel, uniformlyspaced and the interface between fluids is horizontal. The flow here is essentially1D and compares well with the 1D model results, shown earlier in §3.1.2, at leastfor the earlier part of the run. Note that the density and frictional pressure modelsare identical.48(a)(b)Figure 5.1: Two-dimensional simulation of displacement of Example (i)with no eccentricity. a) From left to right, contours of mud mass frac-tion, density, quality, axial velocity and pressure at final time t = 1991s. b) Snapshots of mud mass fraction Y1. White lines show mass stream-lines.49The difference between the 1D and annular models emerges only later in theflow, as the interface nears the top (see the last panel in Figure 5.1b), when aflow instability emerges close to the interface. The parameters of this exampleare such that there remains a significant viscosity gradient between foam and mudover the ranges of quality experienced, i.e. the foamed slurry is more viscous. Asthe displacement progresses the foam density at the interface decreases below thatof the mud and the mean velocity consequently increases above that of the mud.This switch in mechanical stability occurs just after the interface has reached halfway up the annulus; see Figs. 3.1 & 3.6a. But although mechanically unstable, itappears that a threshold density difference needs to be attained before the instabilityis triggered.This is similar to most density-driven porous media miscible displacement in-stabilities, i.e. the Muskat problem, and of course at its heart we have here a mis-cible displacement upwards along a uniform vertical Hele-Shaw cell. When themechanical stability is lost in such systems the destabilizing buoyancy stress stillneeds to overcome a stabilizing viscous effect. For example, with 2 Newtonian flu-ids in a vertical Hele-Shaw, thin penetrating fingers grow ahead of the front when:6Wˆ0(µˆ1− µˆ2)+(ρˆ1− ρˆ2)gˆHˆ2 > 0, (5.1)where fluid 1 is displaced upwards by fluid 2 at speed Wˆ0 and 2Hˆ is the gap width.This makes the threshold explicit. Thus, we expect a similar threshold behaviorhere, although more complex due to the rheology.Figure 5.1a, shows the other flow variables at the final time of Figure 5.1b,we see that the mass fractions, density, quality and velocity have evident pertur-bations. The pressure field does not vary with φ . We should acknowledge thatour imposition of symmetry conditions at the wide and narrow side of the annulus,although convenient for reducing computational time, restricts the azimuthal prop-agation of this instability unphysically. However, the axial propagation is of moreinterest. Behind the unstable region we see that the underlying density gradients inthe foam are stabilizing, so there is no driving mechanism to propagate upstreamof the front.Ahead of the front the mud density is constant and it appears that fingers should50be able to propagate. However, on close inspection it appears that there is a layerof dense mixed fluid that results at the interface, denser than the mud downstream.This mixed layer removes the physical mechanism for fingers to propagate down-stream ahead of the front, hence limiting the unstable region as shown. The originof this higher density mixed region lies with the interpolation of the density. Aswe track advection of both Yg and the individual fluid constituents Y1, Y2, etc., itfollows that the each of these mass fractions must changes close to the displace-ment front. Ahead of the displacement front we have no gas fraction: hence Yg→ 0across the front. Also we have Y2→ 0 across the front. The combination of theseleads to a small band of concentrations over which the mixture density exceedsthat of the mud, i.e. because the slurry density is increasing as Yg → 0. It is hardto see how to rectify this easily, and perhaps also the notion that Y2→ 0 with fixedgas fraction Yg, may not be correct when the two fluids can mix. We leave this fora future model refinement, as the main point here is the onset of instability. Afterinstability occurs, a number of assumptions anyway break down in all Hele-Shawtype displacement models. We shall also see that other effects may also dominateinstability.In Figs. 5.2 - 5.4 we increase the eccentricity of the annulus in the simulationof Figure 5.1, to e = 0.1, e = 0.3 and e = 0.5, respectively. We see an interestingeffect as the displacement progresses. First, low down in the annulus where theslurry density is high we see a steady displacement front. Even with significanteccentricity a steady front may result for sufficiently positive density and viscositygradients between fluids; see [32, 34]. This appears true initially. but by the 3rd and4th time frame of Figure 5.2b we begin to see that the density difference is now nolonger sufficient to displace the mud on the narrow side of the annulus: the frontbegins to lag behind. As the front moves still higher and the slurry density fallsbelow the mud density, we see that the asymmetry of the flow (between wide andnarrow sides) is enhanced by the buoyancy driven instability. Relative to the meanflow, the narrow side fluids move backwards down the well and the front elongatesalong the narrow side. We also observe the onset of density driven fingering in theazimuthal direction, i.e. from the narrow to wide side, as the tail of the displacementextends.51(a)(b)Figure 5.2: Same as Figure 5.1, except the well has eccentricity of e = 0.1.Wide and narrow sides are denoted with W and N.These same effects occur earlier and are more severe, for e = 0.3 and e =0.5. Deep in the annulus the narrow side is not completely displaced and a severe52channel of mud is left behind along the narrow side as the foamed slurry expandsup the wide side. There are more azimuthal instabilities behind the front and thesepenetrate further into the wide side, nearly cutting off the stream of low densityslurry.53(a)(b)Figure 5.3: Same as Figure 5.1, except the well has eccentricity of e = 0.3.Wide and narrow sides are denoted with W and N.54(a)(b)Figure 5.4: Same as Figure 5.1, except the well has eccentricity of e = 0.5.Wide and narrow sides are denoted with W and N.In order clarify more clearly that these instabilities result primarily from thechanges in foam cement density, we have modified this example. Considering55again the example (i) parameters, we now select some representative properties ofthe cement slurry before the interface becomes mechanically unstable. Here wetake the density and rheology values of the foamed slurry when the interface is atξˆ = 150 m and tˆ = 1004 s (as marked by the dashed red line in the fourth panelin Figure 5.1b). Using these values we run a simulation in which the displacingcement slurry is no longer compressible, i.e. Yg = 0 for all time, but with thesesame properties: ρˆ = 1524 kg/m3, n = 1, κˆ = 0.055 Pa.s and τˆY = 7.94 Pa. Figure5.5 shows snapshots of the displacement for four values of eccentricity e = 0, 0.1and 0.5. For e= 0 & 0.1 the displacement front advances at the steady mean speedalong the annulus, displacing perfectly and with no sign of instability. At e = 0.5we also see a steady displacement, but with a very narrow residual channel leftbehind on the narrow side of the annulus. Again there is no instability.56(a)(b)(c)Figure 5.5: Incompressible displacement with cement properties read at ξˆ =150 m and tˆ = 1004 s in Figure5.1. a) e = 0; b) e = 0.1; and c) e = 0.5.57Figs. 5.6 & 5.7 present 2D simulations of displacement in examples (ii) & (iii).(a)(b)Figure 5.6: Two-dimensional simulation of displacement of Example (ii)with no eccentricity.58As discussed in 3.1.2, in example (ii) the entire operation is pumped at a higherpressure, compared to example (i). Hence, the effects of pressure variations aresmaller, frictional pressure effects are smaller and the variation in quality, foamdensity and velocity are smaller. Recall that the inflow quality is fixed at the top ofthe casing (inside), but is now at higher pressure. As a consequence, the foamedcement enters the bottom of the annulus with a higher quality (smaller density)and reaches the density of the mud much lower down in the annulus, compared toFig. 5.1, This is in good agreement with our 1D simulation results in Figs. 3.2& 3.6a. Later, after the density difference between the fluids passes a criticalthreshold we see that the instablity starts to grow. As the foam advances upwards inthe annulus its mean velocity (W0) increases (stabilizing) but the density differencealso increases (destabilizing); see e.g. (5.1). Note also that the total mass flow rateat the surface has changed in example (ii) compared to example (i) (different pˆa).Therefore, the mass streamlines are wider spaced in Figure 5.6.59(a)(b)Figure 5.7: Two-dimensional simulation of displacement of Example (iii)with no eccentricity.Similar to our 1D simulation results, comparing Figs. 5.1 & 5.7, as the in-flow quality is increased from 40% to 60%, the foamed cement enters the annular60section with a density which is already smaller than the density of the mud (seeFigure 3.6b) and this leads to a buoyancy driven unstable displacement happenthroughout the annulus, but starting much earlier. Mass streamlines are againspaced differently compared to previous results because of different total massflowrate. It is notable that the positions at which instability starts in the last fewimages of Fig. 5.6 & 5.7 remains constant and the displacements are stable belowthis heights.To close this section, we now consider a variation of example (i) in which Mode1B conditions are implemented instead of Mode 1A, holding a back pressure on theannulus. This has the effect of pressurizing the entire flow path. Thus, the initialinflow gas fraction is also pressurized, which means that the slurry is compressesless as it travels around the flow path, compared to example (i). Overall the result isthat the compressed foam density at bottom hole is reduced and we retain a lighterweight slurry throughout the annulus. We consider two cases in which the backpressure at the top of the annulus is held at pˆa = 20 atm and 40 atm, respectively.The simulations are slightly more complex numerically as at each timestep we needto iteratively determine the bottom hole pressure in the annulus, so as to attain thetarget pˆa.Figure 5.8 shows colourmaps of the physical variables at the latest time: tˆ =1991 s, as well as displacement snapshots of the mud mass fraction and mass flowstreamlines. The annuli are concentric but here after the density instability startswe see asymmetric channeling along the annulus as the flow develops. This effectstarts deeper in the well when the back pressure is 40 atm and the foam densitiesare lower (Figure 5.9). There remains a threshold for instability, as before, i.e. me-chanical instability is necessary but not sufficient. However, we can see from thesesimulations that the threshold is relatively low for the parameters in these simula-tions, e.g. 50−100kg/m3 is sufficient to trigger these instabilities.61(a)(b)Figure 5.8: Same as Figure 3.1, except a back pressure of pˆa = 20 atm isenforced on top of annulus.62(a)(b)Figure 5.9: Same as Figure 5.8, except a back pressure of pˆa = 40 atm isenforced on top of annulus.63Chapter 6Conclusion6.1 Summary of Foamed CementingThis work outlines development of a mathematical and computational model formiscible displacement of foamed fluids during the primary cementing process. Thebasis of the model is a miscible mixture model in which the foamed slurry is mod-elled as a homogeneous compressible fluid, i.e. the foamed slurry is defined by thelocal mass fraction of gas there is no slip between gas and liquid cement phases.First we have derived a simple 1D hydraulics-based model for the displacementflow. We have also simplified this model to a purely hydrostatic model and usedthe comparison between the two models to show that in many cases the frictionalpressure losses can be neglected, relative to the hydrostatic pressure.Although negligible compared to the static pressures, frictional pressures areimportant in determining the annular flow, and a model for annular displacementflow is one of the main novel contributions of this work. The model derived iscoupled to the 1D model, which provides input into the annulus. Somewhat sur-prisingly, the form of the pressure (or streamfunction) equation is near-identicalwith that for incompressible displacements [7, 27]. The main difference is that,due to the compressibility, we use a mass streamfunction ΨM, instead of the usualvolumetric streamfunction.Numerical solution of the annular displacement model shows that a buoyancydriven instability is triggered at some threshold value of density difference, when64the stabilizing viscosity ratio is also overcome. More clearly, as the foamed slurrydisplaces upwards in the annulus the density falls, eventually dropping below thatof the drilling mud. At some distance above this, the instability starts. Althoughwe have not performed a stability analysis for these two fluids, this threshold typeof instability is reminiscent of miscible porous media displacement instabilities.The instability is triggered lower in the annulus for a more eccentric annulus,as the density difference is progressively less effective on the narrow side of theannulus. Thus, we observe a significant residual mud layer/channel emerge on thenarrow side of the annulus. By comparison with an incompressible displacement,we are able to identify that the root cause of the instability is the decreasing densityof the foamed slurry. Reduction of the effect by application of a back pressure onthe annulus in fact had the effect of reducing the density variation in the foamedslurry and hence triggering the density-driven instability lower in the well.The occurrence of this type of flow instability raises serious questions about theusage of foamed cements. A key motivation for using foamed cements is controlof the hydrostatic gradient in the annulus, i.e. low density. If it is necessary to keepthe slurry density above the mud density to avoid instabilities, this limits the rangeof density reduction. Although there are other benefits to using foamed cements(as discussed earlier), we should also note that there are other lightweight cementalternatives and other ways to achieve similar benefits to foam with use of differentadditives.Specifically, foamed cement appears problematic from the point of view of hy-drodynamic instability. These density driven instabilities are amplified as the foamexpands along the annulus, which self-reinforces the destabilizing mechanism. Al-though here we have used a cement slurry that is more viscous than the mud, this isnot always the case and the uncertainty of foamed cement rheology with pressure(and temperature) means that a viscosity gradient threshold might not be reliable.Even if the foamed cementing operation is undertaken under a controlled pressure,loss of control in this way presents dangers. As in well control, if the annulus isclosed/restricted as the foamed slurry rises, the annular pressure will increase. Ifthe outflow is not restricted then the destabilizing expansion will continue, increas-ing α and eventually leading to coalescence of the gas.656.2 Limitations of the Model & Future DirectionsFinally, we should admit the limitations of our model, which is quite preliminary.These are some future directions for our work:• The impact of temperature variation over the wellbore has been disregardedin our calculations. While, based on (2.3) it can affect the density of thegas. Temperature gradient within the upper layers of the earth (about 100km) is relatively small (15-30◦C/km). As the foam flows down inside thewell, pressure growth increases the density of the gas phase of the foam,while temperature rise causes a drop in the gas density and compensates theeffect of pressure on the gas density, to some extent. It should be noted thattemperature growth decreases the density of the liquids which are assumedto be constant in our model. Including temperature effects and modifyingthe gas density closure are straightforward improvements.• Developing models in order to predict the rheological behavior of foamsis relatively complicated. Moreover, as foams are not always stable undershear, rotational viscometers should only be used with caution for measur-ing the rheological properties of foams. A lot of effort should be made onimproving the cement rheology models, to account for a full range of α andto consider e.g. slip and depletion effects at larger α .• There is a wider theoretical background to develop, underlying pressure-dependent field equations of type (4.37) which arise from both compressiblefluids and α-dependent rheology. Although to date we have worked withthe stream-function formulation, here and previously [10, 27, 33] the (dual)pressure problem might be a better option for compressible flows. Note toothat physically we have other cementing scenarios where compressibilitymay be important (e.g. deep wells with oil-based drilling muds).• For calculating the density of the gas we have used the ideal gas law. Asit was mentioned, the gas which is commonly used for creating foamed ce-ment, is nitrogen (diatomic molecule). Ideal gas law assumptions make it anaccurate equation for monoatomic gases at low pressures and high temper-66atures. Consequently, more realistic equations of state, such as the Van DerWaals equation could be used for calculating the density of nitrogen [38].67Bibliography[1] Canada’s energy consumption and population. https://candaianation.wordpress.com/2014/11/09/canadas-energy-consumption-and-population/.[Online; accessed 1-August-2017]. → pages 1, 2[2] R. M. Ahmed, N. E. Takach, U. M. Khan, S. Taoutaou, S. James, A. Saasen,and R. Godøy. Rheology of foamed cement. Cement and ConcreteResearch, 39(4):353–361, 2009. ISSN 0008-8846. → pages 6, 13[3] P. E. Aranha, C. Miranda, W. Cardoso, G. Campos, A. L. Martins, F. C.Gomes, S. B. de Araujo, and M. Carvalho. A comprehensive theoretical andexperimental study on fluid displacement for oilwell-cementing operations.Society of Petroleum Engineers, paper number 150276, 27(04):596–603,2012. → pages 21[4] O. G. Benge, J. R. McDermott, J. C. Langlinais, and J. E. Griffith. Foamedcement job successful in deep HTHP offshore well. Oil and Gas Journal, 94(11), 1996. → pages 6[5] A. H. Beyer, R. S. Millhone, and R. W. Foote. Flow behavior of foam as awell circulating fluid. Society of Petroleum Engineers, paper number 3986,1972. → pages 17[6] S. H. Bittleston and D. J. Guillot. Mud removal: research improvestraditional cementing guidelines. Oilfield Review, 3(2):44–54, 1991. →pages 5[7] S. H. Bittleston, J. Ferguson, and I. A. Frigaard. Mud removal and cementplacement during primary cementing of an oil well–laminar non-Newtoniandisplacements in an eccentric annular Hele-Shaw cell. Journal ofEngineering Mathematics, 43(2):229–253, 2002. ISSN 1573-2703. → pages3, 4, 20, 21, 36, 43, 44, 45, 6468[8] R. E. Blauer, B. J. Mitchell, and C. A. Kohlhaas. Determination of laminar,turbulent, and transitional foam flow losses in pipes. Society of PetroleumEngineers, paper number 4885, 1974. → pages 17[9] M. Bogaerts, C. Azwar, M. Bellabarba, M. Dooply, and A. Salehpour.Wellbore cementing: an integral part of well integrity. In OffshoreTechnology Conference, 2015. → pages 11[10] M. Carrasco-Teja, I. A. Frigaard, B. R. Seymour, and S. Storey. Viscoplasticfluid displacements in horizontal narrow eccentric annuli: stratification andtravelling wave solutions. Journal of Fluid Mechanics, 605:293–327, 2008.→ pages 21, 66[11] Z. Chen, M. Duan, S. Z. Miska, M. Yu, R. M. Ahmed, and J. Hallman.Hydraulic predictions for polymer-thickened foam flow in horizontal anddirectional wells. Society of Petroleum Engineers, paper number 105583,24:40–49, 2009. → pages 17, 19[12] R. Crook, S. Moore, J. Foreman, and M. Miller. Fully engineeredfoam-cementing process improves zonal isolation. Haliburton EnergyServices Inc. , Drilling contractor, 2000. → pages 8[13] E. P. de Pina and M. Carvalho. Three-dimensional flow of a Newtonianliquid through an annular space with axially varying eccentricity. Journal ofFluids Engineering, 128(2):223–231, 2006. → pages 20[14] L. Ducloue´, O. Pitois, J. Goyon, X. Chateau, and G. Ovarlez. Rheologicalbehaviour of suspensions of bubbles in yield stress fluids. Journal ofNon-Newtonian Fluid Mechanics, 215:31–39, 2015. → pages 11, 12, 14, 15[15] M. B. Dusseault, R. E. Jackson, and D. Macdonald. Towards a road map formitigating the rates and occurrences of long-term wellbore leakage.University of Waterloo, 2014. → pages 2[16] A. Edrisi. Experimental and modeling study of foam flow in pipes with twofoam-flow regimes. PhD thesis, Louisiana State University, 2013. → pages13[17] N. El Khouja, N. Roquet, and B. Cazacliu. Analysis of a regularizedBingham model with pressure-dependent yield stress. Journal ofMathematical Fluid Mechanics, 17(4):723–739, 2015. → pages 4769[18] R. N. Gajbhiye and S. I. Kam. Characterization of foam flow in horizontalpipes by using two-flow-regime concept. Chemical Engineering Science, 66:15361549, 2011. → pages 13[19] B. S. Gardiner, B. Z. Dlugogorski, and G. J. Jameson. Rheology offire-fighting foams. Fire Safety Journal, 31(1):61–75, 1998. → pages 18[20] F. C. Gomes and M. S. Carvalho. Lubrication model for the flow of severalliquids through an annular space with varying eccentricity. Proceedings ofCOBEM, 2009. → pages 20, 21, 22[21] F. C. Gomes and M. S. Carvalho. An asymptotic cylindrical model for wellcementing process. 13th Brazilian Congress of Thermal Sciences andEngineering, Uberlaˆndia, 2010. → pages 20, 21, 22[22] D. J. Guillot, J. Desroches, and I. A. Frigaard. Are preflushes reallycontributing to mud displacement during primary cementing? Society ofPetroleum Engineers, paper number 105903, 2007. → pages 4[23] R. W. Hanks. The axial laminar flow of yield-pseudoplastic fluids in aconcentric annulus. Industrial & Engineering Chemistry Process Designand Development, 18(3):488–493, 1979. → pages 20[24] K. Kopp, S. Reed, J. Foreman, B. Carty, and J. Griffith. Foamed cement vs.conventional cement for zonal isolation-case histories. Society of PetroleumEngineers, paper number 62895, 2000. → pages 5, 6, 10[25] E. W. Llewellin, H. M. Mader, and S. D. R. Wilson. The rheology of abubbly liquid. Proceedings of the Royal Society of London A: Mathematical,Physical and Engineering Sciences, 458(2020):987–1016, 2002. ISSN1364-5021. → pages 11[26] A. Maleki and I. A. Frigaard. Axial dispersion in weakly turbulent flows ofyield stress fluids. Journal of Non-Newtonian Fluid Mechanics, 235(Supplement C):1 – 19, 2016. ISSN 0377-0257. → pages 20, 26, 27, 46[27] A. Maleki and I. A. Frigaard. Primary cementing of oil and gas wells inturbulent and mixed regimes. Journal of Engineering Mathematics, pages1–30, 2017. → pages 21, 36, 41, 43, 44, 45, 46, 47, 64, 66[28] A. L. Martins and W. Campos. A model to design the foam cement job.Society of Petroleum Engineers, paper number 23644, 2(01):43–48, 1994.→ pages 15, 3070[29] R. H. McLean, C. W. Manry, and W. W. Whitaker. Displacement mechanicsin primary cementing. Society of Petroleum Engineers, paper number 1488,19(02):251–260, 1967. → pages 5[30] E. B. Nelson and D. J. Guillot. Intermixing of cementing fluids:understanding mud displacement and cement placement. SchlumbergerEducational Services, 2006. → pages 3, 20[31] M. E. Ozbayoglu, E. Kuru, S. Miska, and N. Takach. A comparative studyof hydraulic models for foam drilling. Journal of Canadian PetroleumTechnology, 41(06), 2002. → pages 17, 18, 19[32] S. Pelipenko and I. A. Frigaard. On steady state displacements in primarycementing of an oil well. Journal of Engineering Mathematics., 46:1–26,2004a. → pages 21, 51[33] S. Pelipenko and I. A. Frigaard. Two-dimensional computational simulationof eccentric annular cementing displacements. IMA Journal of AppliedMathematics, 69:557–583, 2004b. → pages 21, 66[34] S. Pelipenko and I. A. Frigaard. Visco-plastic fluid displacements innear-vertical narrow eccentric annuli: prediction of travelling-wave solutionsand interfacial instability. Journal of Fluid Mechanics, 520:343–377, 2004c.→ pages 51[35] V. G. Reidenbach, P. C. Harris, Y. N. Lee, and D. L. Lord. Rheological studyof foam fracturing fluids using Nitrogen and Carbon-dioxide. Society ofPetroleum Engineers, paper number 12026, 1(01):31–41, 1986. → pages 18[36] A. C. Rust and M. Manga. Bubble shapes and orientations in low re simpleshear flow. Journal of Colloid and Interface Science, 249(2):476480, 2002.ISSN 0021-9797. → pages 11[37] V. Sanghani and C. Ikoku. Rheology of foam and its implications in drillingand cleanout operations. Journal of Energy Resources Technology, 105(3):362–371, 1983. → pages 17[38] R. Span, E. W. Lemmon, R. T. Jacobsen, W. Wagner, and A. Yokozeki. Areference equation of state for the thermodynamic properties of Nitrogen fortemperatures from 63.151 to 1000 K and pressures to 2200 MPa. Journal ofPhysical and Chemical Reference Data, 29(6):1361–1433, 2000. → pages6771[39] P. Szabo and O. Hassager. Flow of viscoplastic fluids in eccentric annulargeometries. Journal of Non-Newtonian Fluid Mechanics, 45(2):149–169,1992. → pages 20[40] A. Tehrani, J. Ferguson, and S. H. Bittleston. Laminar displacement inannuli: a combined experimental and theoretical study. Society of PetroleumEngineers, paper number 24569, 1992. → pages 20[41] P. Valko and M. J. Economides. Volume equalized constitutive equations forfoamed polymer solutions. Journal of Rheology, 36(6):1033–1055, 1992. →pages 18[42] I. C. Walton and S. H. Bittleston. The axial flow of a Bingham plastic in anarrow eccentric annulus. Journal of Fluid Mechanics, 222:39–60, 1991. →pages 20[43] S. T. Zalesak. The Design of Flux-Corrected Transport (FCT) Algorithmsfor Structured Grids, pages 23–65. Springer Netherlands, Dordrecht, 2012.ISBN 978-94-007-4038-9. → pages 27, 4772Appendix ANumerical SimulationParametersFor the two-dimensional results presented in Chapter 5, we have considered an 8in diameter casing inside a uniform 10 in wellbore. The wellbore is assumed to bevertical and 500 m long and it is initially full of drilling mud with the followingdensity and rheological properties: ρˆ1 = 1400 kg/m3, κˆ1 = 0.02 Pa.sn1 , n1 = 1,τˆY,1 = 4 Pa.The mud is displaced with foamed cement, and the properties of the pure ce-ment slurry are as follows: ρˆ2 = 1800 kg/m3, κˆ2 = 0.04 Pa.sn2 , n2 = 1, τˆY,2 = 8Pa.In solving our equations numerically, our main priorities are speed and robust-ness. Hence, for runs of foam cementing displacement simulations we have used100 mesh cells in axial direction and 30 mesh cells in azimuthal direction. Themesh is assumed to be uniform.At each time step, equation (4.36) is solved using slice model to get the massstream function. Using the stream function the velocity field could be constructedfrom (4.26). The fluid mass concentrations are then calculated from (4.27) usingFCT scheme.73