Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Displacement flow of non-Newtonian fluids in an eccentric annulus Pelipenko, Sviatoslav 2002

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-ubc_2002-0200.pdf [ 3.47MB ]
Metadata
JSON: 831-1.0080071.json
JSON-LD: 831-1.0080071-ld.json
RDF/XML (Pretty): 831-1.0080071-rdf.xml
RDF/JSON: 831-1.0080071-rdf.json
Turtle: 831-1.0080071-turtle.txt
N-Triples: 831-1.0080071-rdf-ntriples.txt
Original Record: 831-1.0080071-source.json
Full Text
831-1.0080071-fulltext.txt
Citation
831-1.0080071.ris

Full Text

D I S P L A C E M E N T F L O W OF N O N - N E W T O N I A N FLUIDS IN A N E C C E N T R I C A N N U L U S by SVIATOSLAV P E L I P E N K O MMath (Mathematics) University of Oxford, 2000 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C O L U M B I A January 2002 © Sviatoslav Pelipenko, 2002 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Abstract In this thesis we derive, analyze and devise a method of numerical solution for a Hele-Shaw model of displacement flow of non-Newtonian fluids in an eccentric annulus. The physical problem stems from an industrial process of oil well cementing during the well's construction and successful mathematical modelling and solution of the problem allows for optimization of the process resulting in economic and environmental benefits. Here we outline derivation of the model based on using the long-thin geometry of the annular domain to reduce the flow equations to two spatial dimensions together with the Hershel-Bukley constitutive equations and time evolution equation. Therefrom we obtain analytical solution for the form of the interface in cases of concentric annulus and small annular eccentricity. We proceed to put the problem into its variational and minimization formulation from which we show the existence and uniqueness of a weak solution to the original model. We apply an iterative augmented Lagrangian method to obtain this solution together with a Flux Corrected Transport method for time evolution to arrive at a fully 2-D numerical simulation of the flow. We derive another model using ideas of thin-film and lubrication flows which allows for a quicker prediction of the displacement flow type. We compare the predictions for the flow type based on the lubrication model to those obtained using the 2-D flow simulations. We conclude with a discussion on the significance of results achieved in this work together with relative merits and limitations of the derived models and solutions. ii Table of Contents Abstract ii Table of Contents iii List of Figures iv Acknowledgement vi Chapter 1. Introduction 1 1.1 Process description 1 1.2 Existing literature-review : 2 1.3 Outline of the thesis 4 Chapter 2. Modelling of the flow 6 2.1 Non-dimensionalisation 8 2.2 Model Derivation 11 2.3 Preliminary Results 15 2.3.1 Properties of x : 15 2.3.2 Solution for concentric annulus 17 2.3.3 Perturbation Solution For Small Eccentricity 18 Chapter 3. Analytical results 24 3.1 Variational Formulation • 25 3.1.1 Immiscible fluids with no surface tension 26 3.1.2 Minimization Problem 28 3.1.3 Existence & Uniqueness Results 29 3.2 Solution by augmented Lagrangian method 33 3.2.1 Principle of the method 34 3.2.2 Application of the Algorithm 36 3.2.3 Summary and convergence 38 Chapter 4. Numerical Results 39 4.1 Discretization and method of solution 39 4.2 Convergence 42 4.3 Time advance scheme 47 4.4 Simulation results 49 Chapter 5. Lubrication model 53 5.1 Model Derivation 54 5.1.1 Lubrication assumptions 55 5.1.2 Scalings and approximation derivation 56 5.2 Kinematic equation and interface stability 57 5.2.1 Determining the displacement type 59 iii Table of Contents 5.2.2 Numerical results ' 60 5.2.3 Comparison with flow simulations 64 C h a p t e r 6. C o n c l u s i o n a n d D i s c u s s i o n 7 0 . B i b l i o g r a p h y 73 iv List of Figures 2.1 Schematic picture and a cross-section of the well 7 2.2 Plots of the interface perturbation function g\ for varying density of the displaced fluid P2- Least elongated to most elongated: p2 = 0.5, p2 = 0.6, pi — 0.7, p2 = 0.8, p2 — 0.9, P2 = 1.0, P2 = 1.1- Other parameters fixed at: k\ = 1.0; &2 = 0.9; m\ = m 2 = 1.1; T y I = 1.0; 7V2 = 0.9; pi = 1; 8 = 0; L = 10 22 2.3 Plots of the interface perturbation function g\ for varying yield stress of the displaced fluid TY2- Least elongated to most elongated: ry2 = 0.6, 7y2 = 0.7, TY2 = 0.8, 1Y2 = 0.9, TY2 = 1.0, 7Y2 = l - i - Other parameters fixed at: ki = 1.0; = 0.9; m i = rri2 = 1.1; 7vi = 1.0; pi = 1; P2 = 0.9; /? = 0; L = 10 22 2.4 Plots of the interface perturbation function g\ for varying consistency of the displaced fluid &2- Least elongated to most elongated: k2 = 0.6, &2 = 0.7, &2 = 0.8, &2 = 0.9, k2 — 1.0, fc2 = 1.1. Other parameters fixed at: k\ = 0.9; m i = 7712 = 1.1; T y i = 1.0; TVs = 0.9; px = 1; p2 = 0.9; /? = 0; L = 10 23 2.5 Contour plots of the perturbed stream function $i(c/>, z) in a moving frame of reference. Parameters fixed at: k\ — 0.8; k2 = 1; mi = m,2 — 1.1; Tyi = 1.0; r y 2 = 0.9; p\ = 1; p 2 = 0.9; ^ = 0; L = 10 23 4.1 Schematic picture of the control volume CVij 41 4.2 Plot of Logw of the L2 (solid line) and L°° (broken line) norms of Ap% against the number of iterations together with surface plots of the </>-component of ApJJ after 100 iterations (top to bottom): a) horizontal interface with e = 0.1, b) horizontal interface with e = 0.6, c) slanted interface with e = 0.8 and high yield stress of the displaced fluid. The mesh size here as well as in all of the subsequent examples was taken to be 20 x 80 43 4.3 Plot of Logio of the L2 (solid line) and L°° (broken line) norms of A.PU against the number of iterations together with surface plots of the (^-component of APU after 100 iterations (top to bottom): a) horizontal interface with e = 0.1, b) horizontal interface with e = 0.6, c) slanted interface with e = 0.8 and high yield stress of the displaced fluid. 45 4.4 Plot of Logic, of the L2 (solid line) and L°° (broken line) norms of the filed equations residual Res = V a • S% + fh against the number of iterations together with surface plots after 100 iterations (top to bottom): a) horizontal interface with e = 0.1, b) horizontal interface with e = 0.6, c) slanted interface with e = 0.8 and high yield stress of the displaced fluid 46 4.5 Plot of interface in a concentric annulus at t = 30 together with interface track at 4> = 0 (top to bottom): a) vertical well 8 = 0 (e = 0; pi = 1; P2 = 0.6; Tyi — 1; 1Y2 = 0.8; ki = 1; fc2 = 0.8), b) horizontal well 8 — n/2 (e = 0; pi = 1; P2 = 0.6; ry\ = 1; TY2 — 0.8; k\ = 1; k2 — 0.8), c) inclined well with angle of inclination 8 = 57r/12 and identical rheologies (e = 0; pi = 1; P2 = 0.8; Tyi = 1; TK2 = 1; h = 1; k2 = 1). Inverse power index mi = m2 = 1.0 throughout 51 List of Figures 4.6 Plot of the interface position given by the flow simulation (solid line) and perturbation solution (broken line) together with the track of the interface position on the wide side of the annulus (top to bottom): a) e = 0.1; p i = 1; p 2 = 0.8; Tyi = 1; T y 2 = 0.9; k\ = 1; fc2 = 0.9, b) e = 0.2; p i = 1; p 2 = 0.8; TYI = 1; 7 y 2 = 0.8; ki = 1; k2 = 0.8. Inverse power index mi = m 2 = 1.1 and angle on inclination /? = 0 throughout 52 5.1 Plots of the stream function q and its derivative with respect to <E> versus the interface position (p for different situations. Top to bottom: a) mud channel on the narrow side (e = 0.8; pi = p 2 = 1; T y i = 1; T y 2 = 1-8; k\ = 0.8; k2 = 1); b) steady displacement, no shocks (e = 0.1; p\ = 1; p2 = 0.8; T y i = 1; T y 2 = 0.8; fci = 1; k2 = 0.8); c) unsteady displacement, no shocks (e = 0.1; p i = 0.8; p 2 = 1; T y i = 0.8; T y 2 = 1; k\ = 0.8; k2 = 1). Power index m\ = m 2 = 1.0, angle of inclination (3 — 0 throughout 61 5.2 Plots of the stream function q and its derivative with respect to <& versus the interface position <j> for different situations. Top to bottom: a) shock on the wide side (e = 0.6; Pi = 1; Pi = 0.8; T y i = 1; T y 2 = 0.8; k\ = 1; fc2 = 0.8); b) shocks on wide and narrow sides (e = 0.1; p i = 0.9; p 2 = 1; T y i = 1.2; T y 2 = 0.8; k\ = 1.2; &2 = 0.8). Power index mi — m 2 = 1.0, angle of inclination (3 = 0 throughout 62 5.3 Plots of the stream function q and its derivative with respect to $ versus the interface position 4> for a slumping flow situation (e = 0.6; p i = 1.2; p 2 = 0.6; T y i = 1.2; T y 2 = 0.6; k\ = 1.2; k2 = 0.6). Power index m\ = m 2 = 1.0, angle of inclination (3 = 0 62 5.4 Contour and relief plots of the differential velocity WitW — Wi:n for varying density p 2 and yield stress T y 2 of the displaced fluid and other parameters fixed (given in text above). On the left plot, darker shaded region corresponds to predicted unsteady flow and lighter to predicted steady flow 63 5.5 Contour and relief plots of the differential velocity Wj ) l u — Wi, n for varying density p 2 and consistency k2 of the displaced fluid and other parameters fixed (given in text above). On the left plot the darkest shade corresponds to predicted unsteady displacement, the lighter shade to predicted steady displacement and the lightest shade to predicted slumping (steady) displacement 64 5.6 Contour and relief plots of the differential velocity Wj j U ) — Wi,n for varying density of displaced fluid p 2 and eccentricity e of the annulus and other parameters fixed (given in text above). On the left plot the black region corresponds to predicted mud channel on the narrow, the lighter shade to predicted unsteady displacement with fluid yielded on the narrow side and the lightest shade to predicted steady displacement. 65 5.7 Top: plot of the differential velocity W^w — W^n calculated using the lubrication model versus the density of the displaced fluid p 2 with other parameters fixed (e — 0.4; k\ = 1.0; k2 = 0.8; m i = m 2 = 1.0; T y i = 1.0; T y 2 = 0.8; p i = 1; (3 = 0). On the plot, solid circles indicate the values of density for which 2-D flow simulation indicated an unsteady displacement and empty circles where a steady displacement was achieved. Below are plots of the interface position at the wide side of the annulus (<p = 0) from the simulated flow in frame of reference moving with the average speed flow for (to to bottom, left to right): a) p 2 = 0.4; b) p 2 = 0.5; c) p 2 = 0:6; d) p 2 = 0.7; e) p 2 = 1.0; f) p 2 = 1.1 . . . . 66 5.8 Right: plot of the interface position at the wide side of the annulus from the simulated flow in frame of reference moving with the average speed of the flow for p 2 = 0.7 with interface evolving from an initially elongated configuration (plotted on the left). Other parameters fixed at the same values as in Fig. 1.7 67 vi List of Figures 5.9 Top: plot of the differential velocity WitW — Wi,n calculated using the lubrication model versus the consistency of the displaced fluid k2 with other parameters fixed (e = 0.3; ki = 1.0; mi = m 2 = 1.0; ryi = 1.0; ry 2 = 0.8; pi = 1; p2 = 0.95; 3 = 0). On the plot, solid circles indicate the values of consistency for which 2-D flow simulation indicated an unsteady displacement and empty circles where a steady displacement was achieved. Below are plots of the interface position at the wide side of the annulus (<j> = 0) from the simulated flow in frame of reference moving with the average speed flow for (left to right): a) k2 = 0.4; b) k2 = 0.6; c) k2 = 1.2 68 5.10 Plots of the velocity at the wide side of the annulus (<fi = 0) as given by the lubrication model (solid line) and from 2-D flow simulations (crosses). Left to right: a) density of displaced fluid p2 between 0.9 and 1.3 and other parameters fixed at e = 0.4; k\ = 1.0; k2 = 0.8; mi = m2 = 1.0; ry\ — 1.0; ry2 — 0.8; p\ = 1; 3 = 0 b) consistency of the displaced fluid k2 between 1.0 and 1.8 with other parameters fixed at e = 0.3; k\ — 1.0; mx = m2 = 1.0; Tyi = 1.0; Ty2 = 0.8; pi = 1; p2 = 0.95; 3 = 0 69 vii A c k n o w l e d g e m e n t This research has been supported by Schlumberger Oilfield Services and by NSERC Canada through C R D project 245434. Financial support of the Pacific Institute of the Mathematical Sciences is also gratefully acknowledged. I would like to thank most of all my supervisor, Dr. Ian Frigaard for his invaluable help and guidance. The readers for this thesis were Dr. Brian Wetton and Dr. Anthony Peirce, whose comments and endorsement are greatly appreciated. v i i i Chapter 1 I n t r o d u c t i o n 1.1 Process description In construction of oil and gas wells it is necessary to cement a series of steel casings into the well as the depth increases. This is done both to support the wellbore and to provide zonal isolation throughout the length of the well, ensuring a hydraulic seal on the outside of the casing. The primary cementing process proceeds as follows. After a new section of the well is drilled, the drill pipe is removed from the wellbore, leaving the drilling mud behind. A section of the steel casing is then inserted into the hole, leaving a gap of ~2 cm between the outside of the casing and the rock, i.e. the annulus. Although centralizers are sometimes fitted to the outside of the casing to prevent the steel tubing from slumping to the lower side of the wellbore, it is very common that this annular gap is eccentred. After the casing is in place, a sequence of fluids is pumped down the inside of the tubing, reaching the bottom of the well and coming back up through the annular gap. Typically, a wash or spacer fluid is pumped first, displacing the drilling mud left over both on the inside of the tubing and outside in the annulus, followed by one or more cement slurries. Drilling mud follows the final cement slurry pumped and circula-tion is stopped with a few meters of cement left at the bottom inside the casing. The cement is then allowed to set. Finally the drilling resumes through the remaining part of cement inside the casing and further into the underlying rock. A successful cement job results in removal of mud and spacer fluid from the annulus by the cement slurry. Unfortunately, mud is sometimes left behind in parts of the annulus and as the 1 Chapter 1. Introduction cement sets, water is removed from it, producing a porous channel along which liquids and gasses from the surrounding rock formations can migrate upwards. This may result in negative safety and environmental consequences as well as loss of productivity. To rectify the problem, expensive squeeze cement jobs are required. A primary cement job design, in the context of this thesis, consists of a sequence of fluids, volumes and flow rates, with specified fluid densities and rheologies. To aid annular displace-ment during cementing, the rheologies and densities of the spacer fluid and the cement slurry can be modified within the constraints of maintaining well security. The pump schedule regulating the rate of flow through the gap can also be varied. 1.2 Ex i s t ing literature review A n overview of primary cementing and a description of common industrial practices used is presented in [18, 41]. Methodology for cementation described there is primarily based on expe-rience and several "rules of thumb" for improving the mud displacement are given, supported by physical considerations. Important feature of the displacement flow such as possible presence of a mud channel on the narrow side of the annulus has been identified and examined in [30], where a hydraulic approach is adopted and in [19, 20] together with field evidence. Techniques derived from field experience and supported by large-scale experiments for detection and re-moval of the mud channels and unyielded cement slurries are outlined in [3, 36, 37]. Hydraulic reasoning has been used in the majority of the industrial literature on the subject, leading to several systems of design rules for a successful cementing job as in [23, 27, 28, 38]. Generally, these rules are applicable to near-vertical wells and state that the flow must be suf-ficiently high to avoid a mud channel on the narrow side of the annulus and there should be a positive gradient of Theological parameters and densities (so that each fluid generates higher factional pressure at the wall and is heavier than its predecessor). One such system of rules currently used in well construction is the W E L L C L E A N system, developed by Schlumberger-2 Chapter 1. Introduction Dowell, [8, 14, 18] with examples of application given in [6, 25]. Whilst these approaches are based on a number of physical truths and have been successfully applied to a variety of cementing jobs, their predictions are generally conservative and they experience problems with predicting flows in highly deviated or horizontal wellbores where new problems arise [9, 24]. More comprehensive and detailed solutions fully modelling the flow in 3-D have been computed recently [2, 44]. It is noted however that computational requirements for simulating 3-D flow over the scale of the wellbore are prohibitive and only an annular cross-section of the domain is usually considered. Here we follow the derivation of the model given in [15] that allows for generality in terms of well geometry and a level of detail that the more restrictive hydraulic design rules do not provide. At the same time it makes sufficient number of justified simplifications to reduce the problem to 2-D, thus greatly reducing the computational load required to model displacements over the scale of the wellbore. The general idea is based on averaging in the radial direction reducing the problem to 2-D, in axial and azimuthal, directions in a manner similar to that applied to computing flows in Hele-Shaw cells. The first appearance of this idea, in fact deriv-ing from a porous media flow analogy is in [29]. The current work is based on [15], which has evolved from [4, 5, 10, 13]. Limited experimental validation of the approach is given in [42, 43]. The method of numerical solution for the resulting flow equations used in [15] requires a few assumptions and simplifications additional to those made in derivation of the model. This is prompted by the requirements of speed and robustness for the flow simulator so that fre-quent, quick and trouble-free re-runs of cementing displacement simulations can be made and a near-optimal cementing job design can be found. Here however, we attempt to solve the flow equations without making any additional assumptions or simplifications to the model beyond putting it into variational formulation. We accomplish this by using the augmented Lagrangian method for solution of variational elliptic equations, outlined in [12, 17]. This method has been applied successfully in a simpler situation of a one Bingham fluid flow in an eccentric annular 3 Chapter 1. Introduction cross-section in [22]. Later in this thesis we derive a simplified model for the flow, allowing for quick determina-tion of the displacement type based on lubrication-type assumptions. Essentially we assume the displacing fluid to be "shooting up" on the wide side creating an elongated interface and resulting in several simplifications to the flow equations. By looking at whether the interface elongates or not we conclude whether or nor the given situation would result in an unstable or stable displacement. This type of approach is broadly similar to examination of fingering effects in a non-Newtonian Hele-Shaw displacement flow, which has been studied recently in [26, 7]. It should be noted however that in the case of an eccentric annulus we. have an additional effect of increasing wall friction towards the narrow side when compared to a concentric annulus, or equivalently to an Hele-Shaw cell with parallel plates. Although not studied extensively in this thesis, the effect of this seems to prevent fingering at the interface in most situations causing the flow instabilities to manifest themselves instead by creating one long "finger" at the wide side of the annulus. This comment also applies when comparing the lubrication model presented here to the earlier works on interface instabilities in Hele-Shaw type geometries such as [31, 39, 40] or to interface instabilities for porous media displacements [21, 32, 33, 34, 35]. Also here we consider mostly "stable" displacements, whereas the Hele-Shaw and porous media flow studies listed focus on unstable displacements. 1.3 Outline of the thesis In this thesis we analyze and obtain solutions for a model of the laminar displacement flow of Herschel-Bukley fluids in an eccentric annulus. The outline of the thesis is as follows. In Chapter 2 we proceed to describe the derivation of the model closely following work previ-ously done in this area. Classical dimensional scaling together with averaging across the annular gap analogous to that done in a Hele-Shaw displacement model, lead to a system of equations in two spatial dimensions. These consist of a coupled system of linear advection equations for the 4 Chapter 1. Introduction fluid concentration and a nonlinear Poisson equation for the stream function. Thereafter, we obtain some preliminary analytic results following directly from the formulation of the model - such as expressions for the steady-state interface for concentric annulus and a perturbation solution for small eccentricity. In Chapter 3 we show existence and uniqueness of the weak solution to the nonlinear stream function equation. We accomplish this by putting the problem in variational formulation, con-verting it to a minimization problem and using standard results of functional analysis. Then we apply the augmented Lagrangian method to obtain a solution to the minimization problem through an iterative algorithm. We simplify the steps of the resulting algorithm and state its convergence properties. In Chapter 4 we detail the discretisation and numerical solution of the equations derived in the previous chapter. We examine convergence of the algorithm and compare the results with analytic and perturbation solutions given in Chapter 2. In Chapter 5 we derive a simplified model for evolution of the interface using methods ap-plied to lubrication and thin-film flows. In this way we effectively reduce the process dynamics to a one-dimensional model. It allows us to give a criterion for the type of fluid displacement occurring without actually simulating the entire flow. We give some results based on this model and compare them with those obtained from the 2-D flow simulation. 5 Chapter 2 M o d e l l i n g o f t h e f l o w A local cylindrical coordinate system (r, 0, £) is used to describe the well geometry; £ measures distance along the central axis of the casing r = 0, (i.e. £ is the measured depth, but measured upwards from the bottom of the well). Wells are typically inclined to the vertical and the inclination angle is denoted /?(£)• The local cross-section of the well, outside the casing, is assumed to be that of an eccentric annulus, with inner radius f j ( £ ) , equal to the outer radius of the casing and outer radius r 0 (£) equal to the inner radius of the hole (or previous casing). At each depth £, the mean radius r a (£) and the mean half-gap width d(£) are denned by: ra(0 = \[ro(i) + m)}, d(£) = \[ro(i)-h(0}- (2.1) As well as inner and outer radii, the displacement e(£), of the two centres of the two cylinders is given. The following three geometrical assumptions are made: (i) that the cylinders do not touch, (e(£) < 2ci(£)); (ii) all variations in the cross-section geometry and inclination, axially along the wellbore, are slow; (iii) the weight of the casing acts in such a way that the narrow side of the annulus will be found on the lower side of the well. The flow is assumed to be laminar and incompressible. The Navier-Stokes equations together 6 Chapter 2. Modelling of the Bow mmtmMm&szmmmmA Figure 2.1: Schematic picture and a cross-section of the well. with the incompressiblilty condition in cylindrical coordinates are P lot j I d , . 1 d . d 1 d... . 1 d . d . f e e c9p v = at r 2 c9r 1 (9 1 dp , 1 0 . 0 . dp _ 1 <9 r _ . 1 dv dw rdr1' r d6 Qf (2.2) (2.3) (2.4) (2.5) where u = (u,v,w) is the velocity, p is the pressure and g = (gr,go,gg) is acceleration due to gravity. We now assume that a sequence of K fluids is pumped around the flow path and concentration of each individual fluid component is denoted ck, which is modelled by an advection-diffusion equation: ~o4  + ~r ~dr[™Cfcl +  + ~^Ck^  = V • [Dk(c,u)Vck], (2.6) rd6l dC where 2~2k=i Cfc = 1; C = ( c i , c 2 , . . , c#) a n d the rheologica l parameters together w i t h the dens i ty d e p e n d u p o n the concent ra t ions o f f luids present. 7 Chapter 2. Modelling of the flow 2.1 Non-dimensionalisation Here we non-dimensionalize and scale the equations above with the aim of deriving a two-dimensional model in azimuthal and axial directions. We eventually eliminate the radial de-pendency by averaging across the annular gap. Scalings: Let the axial coordinate of the top of the well be £ = £<z with the bottom of the well at £ = t^bh and thus length of the well in dimensional units Z = t\tz — tbh- Then define 1 f&z , „ i ritz rl = ^ ra(0 d£, 5* = -r I 5(0 d£. (2.7) being the average mean radius and global measure of narrowness of the annulus respectively. Rescale the axial and azimuthal coordinates £ and <fi in a natural way: 7 £ = ^ # , </>=-• (2-8) nr* 7r Also, in each cross-section define local average radius, r a (£ ) , local annulus eccentricity, e(£), and local measure of the narrowness of annulus, <5(£) by: «0 = - % r . ( 0 - ^ , ,(9 = i | . (2.9) ra(0 ra 2d(£) Then, the centreline of the annular gap is at r = ra(£)rc(c/>, £), inner wall at r = n (</>,£) = r a ( O M 0 , O - S(Oh((p,0] and outer wall at r = r*(ci,0 = ra(£)[rc(<M) + <5(OM0,£)], w h e r e for small <5(£) the following expressions define rc(cp,c;) and /i(</>,£): rc(<M) ~ l + e ^ ^ c o s T T ^ + t e ^ ^ ^ s i n T r ^ + O ^ 3 ) , (2.10) /i(0,O ~ l + e(Ocos7r0-[e(O]25(Osin7r</» + [e(Oc5(Osin7r0]2 • +0(<53). (2.11) Assuming sufficient narrowness and-uniformity of the annular gap so <5(£) ~ S* <C 1 we discard terms 0(<5(£)) above, taking: • r c ( 0 , £ ) = l , ^ , O = l + e(Ocos7r0. (2.12) Chapter 2. Modelling of the flow-Now scale the radial coordinate relative to distance to the centerline of the gap, as follows: r-ra(Z)rc(<l>,Z) r - ra(g) , , y = T . = — ^ 7 — I 2 - 1 3 ) Thus the outer and inner walls are given by: v = ± m ( ) = ± m i M ^ h ^ . p.14) To scale the velocities we define a typical cross-section of the annulus A* = 47r<5*[f*]2 and a scale for the flow rate Q* = m.a.-x.^Qpump(i). Then, the axial and azimuthal velocities are scaled with w* = Q*/A* and radial velocities are scaled with W*5*/-K. Finally, define dimensionless flow rate and dimensionless time by t = ^i, = (2.15) irr* . Q* To scale the fluid properties, pressure and deviatoric stress we note that the characteristic scale for the rate of strain is 7 * : 7* = ^ , • (2-16) and we use this to define scales for shear stress, viscosity and pressure, as follows: f* = m F [ f f c , y + « f c t f T * ] . A* = T ; . P*  = Ijr> (2-17) fc <y* 0 where for each fluid k f^y is the yield stress, kk is the consistency, rik) is the power law index and each fluid is assumed to obey Herschel-Bulkley constitutive laws. A l l shear stress components are scaled with f* and dimensionless rheological parameters are defined by: Tk,Y=T*TktY, kk(j*)nk = ft*Kk, nk = nk. (2.18) The fluid densities are scaled with p*\ p* = max[pfc]. (2.19) fc Scaled equations: Applying the scalings above to the momentum equations (2.2 - 2.4) and neglecting terms O(^jr) (where 6*/ir is the ratio of radial to azimuthal length-scales and is 9 Chapter 2. Modelling of the How assumed to be small): dp dy' 1 dp d psin.3 sin ircf) st* _ _dp d _ pcosd U ~ dt + dyT*v St* ' (2.20) (2.21) (2.22) where St* is the Stokes number: st* = ^ (2.23) p*g*[r*6*]2 p*g*f*5*' The leading order scaled momentum equations above describe a bidirectional shear flow through a slot of width 2H((f),£). Now, the leading order expressions for rates of strain are: dw „,5* (2.24) where n is the effective viscosity. Substituting this into the scaled equations above and assuming that the velocity field is approximated by a slot velocity field s = (vs,ws) symmetric about the centreline of the gap we obtain: dvs d_ dy d_ dy dy dws dy _ 1 dp p sin 3 sin 7r</> Tad~4> ST* ' _ dp p cos 3 bX+ St* " (2.25) (2.26) where p = p((j>,€,t) is now independent of the radial coordinate. Scaling the incompressibility condition (2.5) and retaining leading order terms again: du . 1 dvs . dws + = 0. (2.27) dy ra d(f) d£ To eliminate the radial velocity u we average across the gap width and use the no-slip condition at the walls to give: f) Ft (2.28) £[ lW| + | [ r . m | - 0 , where 1 fH 1 rH = #y v*dv> "tt't'Q= itJQ  W s d y >  (2-29) 10 Chapter 2. Modelling of the flow are the slot velocities averaged across the gap. Thus, to satisfy (2.28), we introduce a stream function ^ as in an incompressible 2-D flow: raHw = 9 0 ' u- ^ (2.30) Next, we scale the concentration advection equation (2.6), retain the leading order terms and integrate through the gap as above. Noting that the rescaled diffusion terms appearing on the right hand side are typically of a negligible magnitude for the dimensions of cementing jobs considered, we obtain [Hrack] + [Hv cfc] + [Hraw cfc] = 0, (2.31) where Cfc denotes the gap-averaged fluid concentration. Finally, using the rheology and density scalings above we rescale and reduce the consitutive laws: T<t>y = T]~bS ^ | r ' > T I Y = V dws dy \T\ > Ty(c) 7 = 0 |r| < ry (c ) , (2.32) (2.33) (2.34) dv3 + dws where |r | = [T£, + r 2 / / 2 , 7 = the fluid is yielded, the effective viscosity 77 is: 1/2 is rate of strain and in the regions where 77 = 77(0,7) = K ( C ) 7 ( " ( c ) - 1 ) + TY(C) (2.35) 2.2 Model Derivation By integrating (2.25) & (2.26) three times with respect to y, to extract slot-averaged velocities we get: 1 dp p sin 8 sin 7rc/> w ra d<f> dp p cos 8 St* H r H H rH H Jo ly y ^  dydy, H fo jy v(y) dydy, (2.36) (2.37) 11 Chapter 2. Modelling of the Bow-Note that from (2.20) p is independent of g, and that we have assumed that concentration does not vary across that gap. From this it follows that the vector of averaged velocities (v,w) is parallel in the (</>,£) plane to the vector G: N ( N N N _ ( 1 dp psmBsimrcf) dp pcos8\ G = (G+,Ge)=(--^ + ^ . - - - — - J , (2.38) which is the vector of modified pressure gradients. Hydrau l i c assumption We have from the definition of G, scaled consitutive laws and the reduced scaled momentum equations (2.25) &; (2.26): | - T = - G , (2.39) where r = ( r^ , r^y) satisfies (2.32 - 2.34). This is an equation for a Poiseuille flow between two plates distance 2H apart. By integrating this with respect to y, using the no-slip conditions at the walls and substituting from the constitutive assumptions, we obtain the following relation between the modified pressure gradient G and the rate of flow through the slot H\y2 + I Z J 2 ] 1 / 2 = F(G): ' 0 HG < Ty, (HG-Ty)(\m + l]HG + TY) \HG - Ty F(G) = { K (2.40) HG > Ty, G2~(m + l)(m + 2) where n = 1/m. Now to simplify the notation we introduce the following modified gradient and divergence operators: V o - f - i - ^ ^] V a - 1 • ^ MS)d<t>'dU' ' a * ra(0 dct> ' 6T for arbitrary q and q = (q^, q{). We note that H\s\ = | V a * | , so that (2.40) becomes | V a * | = F(G). Substituting back into (2.36) & (2.37) we have'when | V a * | ^ 0: • _tm I aw ZJL _ G± KSS _ <h ( 2 . 4 1 ) | V a * |  ' | V a * | G From this we can derive a formulation for the pressure field or the stream function. The pressure field is however indeterminate in the unyielded regions where IVa^l = 0, so we proceed with deriving the stream function formulation. 12 Chapter 2. Modelling of the flow Stream function formulation To make explicit the effects of the yield stress on the flow and avoid singularity in equations for the stream function at | V a * | = 0 we define x by: Ty X = G H substituting into (2.41) for G in terms of x-dp rapsinBsinir<f> _ rQ(x(|V a*|) + r y / H ) d V -raGcj, d<t> St* |Va* | dp pcosd x( |V 0 *|) + TY/H 1 d% (2.42) (2.43) (2.44) d£ St* | V a * | rad<l>' valid in the regions where the fluid is yielded. Cross-differentiating to eliminate the pressure dependence we finally obtain: Va-S = - / , (cj),£) £ (0,1) x (0,Z), (2.45) where S = (raG^-raGf) = ra x(lVq*!) + T y / # |Va* | raTY H ' and / is given by: | v a * | = o | S | < rap(c) cos 3 rap(c) sin 3 sin ir(f) rqTy H ' (2.46) (2.47) (2.48) S t * ' St* Finally, substituting for % into the expression for F(G) above we obtain a relation between x and I V Q ^ I if required: |Va* | = < 0 ,m+l [ Km(m + 2) (x + ry/Hf X + (m + 2 ) T y (m+l)H\ X < 0 , x>o. (2.49) Boundary Condi t ions On the wide and narrow side we have: * ( < U , t ) = o, ¥(U,t) = Q(t), (2.50) the value of stream function on wide side is taken to be 0 arbitrarily and the condition on the narrow side arises thus from the incompressibility of the fluids pumped. For conditions on the 13 Chapter 2. Modelling of the flow top and the bottom ends of the annulus, we assume the flow there is in the axial direction only, so d$> — (4>,Z,t)=Q, (0,0,0 = 0. (2.51) The procedure for this concentration-dependent model is to calculate the stream function from the original concentration function by solving the equations above, then extract the slot-averaged velocities by differentiating the stream function. Thereafter we compute the evolution of the concentration from the advection equation (2.31). Immiscible fluids model Instead of the concentration-dependent model above we can also consider the displacement as that of immiscible fluids without surface tension, i.e. since the concentrations are simply advected by (2.31). We show in Chapter 3 that both these models lead to the same variational formulation. Thus, without loss of generality, if we consider a two fluid displacement (with fluid 1 dis-placing and fluid 2 being displaced) and following identical steps as in the above concentration-dependent model we obtain: V A • S i = 0, for x e t t i (2.52) V A • S2 = 0, for x € fi2, (2.53) with Si and S2 defined as in (2.46- 2.47) in each of the fluids and £li,£l2 are the regions where correspondingly fluid 1 and fluid 2 is present. Denote the interface between two fluids by £ = h(<f>, t). Then we require the pressure p and the stream function $ to be continuous across the interface, or b]? = 0; ml = 0, ' (2.54) where denotes the difference in q between fluid 2 and fluid 1 across the interface. The kinematic condition at the interface gives dh v dh _ ,n . dt rad<f> 14 Chapter 2. Modelling of the flow which is the equation for the advection of the interface, replacing (2.31). From the continuity of p and ^ we have that their directional derivative in the direction tangential to the interface t = (1, § | ) is continuous, so: [t • V ap]? = [t • V a * ] ? = 1 dp dpdh ra d<f> <9£ d<j> i_a* av_dh ra d<p d£ d(f> = 0 0. Substituting for components of pressure gradient from (2.43 & 2.44) we obtain: |2 X + ry/H 9* _ ]_d^_dh <3£ ra d(j> d(j> p sin 0 sin 7r</> p cos 0 dh 0 (2.56) (2.57) (2.58) | V a * | J \d£ ad4>d4>) ' \ St* St* d4>, as the pressure continuity condition (2.56) and note that (2.57) is equivalent to the continuity of the normal component of velocity. The latter is required for (2.55) to be valid. 2.3 Preliminary Results It appears that obtaining a solution for the above system of equations for a general case will require use of numerical methods - for example by applying the augmented Lagrangian method to a finite volume approximation for the weak formulation of the problem as we have done in Chapter 3. However, some analytical results for simple cases can be derived and from these, the trends of varying model parameters can be deduced. These results can aid both in validating the numerical solutions for the general case and provide some insight on the qualitative physical level as to the behavior of the flow. 2.3.1 Properties of x The relation between the modified pressure gradient and the rate of flow from (2.49) is given implicitly by X < 0 , 0 Ax m+l m + 2 „ m + l X > 0 , (2.59) where x — G — ry /H with G being the modulus of the modified pressure gradient defined as above and (2.60) rrm+2 A = - £ > 0; Km(m + 2) 15 Chapter 2. Modelling of the flow 1. Asymptotic behavior in the yielded region % > 0: As | V a * | - 0, I V a * | ~ ^ 1 X ~ I V ^ I 1 / ^ 1 ) - D y 771 T J- / As | V a * H oo, | V a * | ~ A x m =^xHV a t f | 1 / r n (2.61) 2. Derivatives with respect to | V a * | for x > 0: dX H V a * 0 | V a * | V 3. Bounds for x - l > 0 (2.62) * A { ^ ) X M > ( 2 - 6 3 ) as (x + )B)/(X + B) is decreasing and bounded by (m + 2)/(m + 1). Thus / m + i \ V™ | v < , , t | 1 " " ( 2 ' 6 4 ) Also, for the upper bound we have: A m+l | V Q * | > • (2.65) X + 5 f Axm+1/(2B) x<B, > I . (2.66) [ A X m / 2 x > B, Rearranging, we obtain: ( ^ i / O + i ) | v ^ | 1 / ( m + 1 ) \va^\<ABm/2, X < < (2.67) I (I)17™ | V a * | 1 / ? n | V 0 * | > ABm/2, 4. Concentric annulus In case of a concentric annulus with eccentricity e = 0 and constant unit flow rate Q(t) = 1 we have X = x(li 1) satisfies: 16 Chapter 2. Modelling of the Bow in the yielded region. Also differentiating (2.59) with respect to | V a ^ | and H we obtain, correspondingly: x'(M) = dx •(1,1) = K T O ( m + 2 ) ( X ( l , l ) + T y ) 3 dx m X ( l , l ) m ( x ( l , I ) 2 + 2^±f X ( l , l)ry + ^ r Y ) (m + 2)(x( l , l ) + T y ) 3 • XH(U) - d H M - T V m ( x ( 1 ) 1 ) 2 + 2 ^ ± 2 x ( 1 ) 1 ) 7 y + ^ T 2 r (2.69) x ' ( i , i ) x m ( i , i ) Km (2.70) thus giving the relation X H ( I , I ) - TV = -which will be used later. 2.3.2 Solution for concentric annulus j Here we obtain a solution for the case of zero eccentricity and constant rate of flow, so e = 0, H = 1 and without loss of generality Q(t) = 1. We consider a problem with only two fluids and adopt the immiscible fluids model with continuity conditions at the interface described in (2.55 - 2.58). First, introduce frame of reference (z,<j>) where z == £ — t is moving with the average speed of the flow. Then, writing z = g((j>,t) = h(</>,t) — t for the position of the interface, the continuity conditions (2.56 - 2.58) become: |2 d$ d§_dg_ d<j) dz d<f> = 0 X + TY ( h D — dz (dcp 'dcp + p sin 3 sin ir(j) p cos 3 dg = 0, (2.71) (2.72) J V a ( * + 0)iy \dz v<90 'd<t>J \ St* St* d<f>, at the interface, where $ = \& — </> is the stream function in the moving coordinates. Now, $ = 0 clearly satisfies the first interface condition above. It also satisfies the field equations as S =  X^ V^ ^'TY V a \ I / is then constant in each fluid. The second interface condition is satisfied if pcos/31 X ( l ) + r y + St* ]2dg psin 3 . 1 0 0 St* sin(7r</>), (2.73) which gives the shape of the interface. So, for a vertical well we have 3 = 0 giving g constant in (f>, so the steady-state interface is horizontal as we would expect. For a horizontal well, 17 Chapter 2. Modelling of the flow integrating the equation above we obtain # ) = - , t o , , : " h l 2cos(7r0) (+const.) (2.74) i w?. *sr [X(l) + Ty\j Note that if the two fluids have identical rheologies in the case of horizontal well we can only obtain a steady state if they also have identical densities, i.e. are the same fluid from the point of view of our model, in which case the position of the interface is essentially undefined. One other case we may consider is that of an inclined well with two fluids with same rheologies but different densities. Then, the steady-state interface is given by g((f,) = -itan/3cos(7rtA) (2.75) 7T Remark In Chapter 3 we show that both the concentration-based model and the immiscible fluids model have the same weak formulation and that formulation has a unique solution. Thus, since we have showed above that $ = 0 satisfies the equations of the immiscible fluids model for a concentric annulus, we can conclude that it is also a unique solution and would satisfy the corresponding concentration formulation, i.e. without any intermediate concentrations. 2.3.3 Perturbation Solution For Small Eccentricity Here we use standard perturbation methods to find the stream functions and shape of the steady-state interface in case of small eccentricity 0 < e << 1 of the annulus. As in the case above of a concentric annulus, we assume a constant flow-rate Q(t) = 1 as well as small constant eccentricity e « 1 and constant mean radius ra — 1 (thus we have V a = V ) . Also, for the sake of simplicity we assume the top and bottom of the well are at z = ±L. We assume that the unperturbed steady-state interface is z = ego((p) in the moving frame of reference with go((f>) = 0(1). Note that from the interface condition (2.73) this requires d<f> [St*X(l) + St*TY + p cos 0\\ W so we need either a nearly-vertical well with the angle of inclination (3 small, or the rheological difference much greater than the density difference between the two fluids, to support the nearly-azimuthal interface. 18 Chapter 2. Modelling of the how-Field Equations From the above we have that the concentric (unperturbed) steady state solution is given by <I> = 0, g = ego(4>), where $ is the stream function in the moving frame and z = g(4>) is the position of the interface in the moving frame of reference. In case of small eccentricity we seek a perturbed solution of the form $(z,c/>) = 0 + e$i(z,</>) + e2$2{z,<t>) + -g(<j)) = eg0{4)) + egi{<f>) + e2 g2{(j>) + ... From this we obtain: V * = V$ + V$* = (1 + eCOS(TT0) + e$i^ + e$i,z + •••), where <£* = </!> + (e/ir) sin7r0. Using the binomial expansion: , |V$ + V$*| ~ 1 + e(cos(-K(f>) + + 0(e2) giving X ( | V $ + V**| , H) + ry/H ~ 1) + e(cos(7r</>) + $i,^)x'(l, 1) +ecos(7r</))xH(l, 1) + TY(1 — ecos(ir(j))) + 0(e2). Thus X + TY/H |V$ +V$* (V$ + V<T) ~ ([x(l, 1) + TV] + e{(cos7T0 + 1) + COS7T(/>XH(1,1) -TycosTT^} , [x(l,l)+ r y ] $ i i 2 ) . Substituting into the field equations, at 0(e) in each fluid we obtain X ' ( l , + (x(l, 1) + 7 Y ) $ M * = (x'(l, 1) + XH(h 1) - Ty)7rsin7f(/> in each of the fluid domains Q\, Q2. (2.77) (2.78) (2.79) (2.80) (2.81) (2.82) (2.83) Boundary and Interface conditions At <f> = 0,1 we have $ i = 0 by the choice of \&*. At z = ± L we have = 0. The continuity conditions at the interface become: l 2 + e$i = 0 (X + 7 V ) ( * 1 , 2 0(0o + S i h , psin/3sin7rc/) pcos/3 d(go + 51) .90 5t* St* <9(/> Ji (2.84) (2.85) 19 Chapter 2. Modelling of the Bow on z = 0. Note that in deriving (2.84) and (2.85) we linearize both about about the same solution and with respect to the unperturbed domains Q i , f22- Note that from the interface advection equation it follows that ^ = 0 only if $i^(c/>, 0) = 0 in both fluids. Thus for steady states we impose $(0,0) = 0. Solut ion of the F i e ld Equations First homogenize the field equation (2.83) by letting U(4>,z) = c h ( < M + - (1 - **(M)-M sinTr^ (2.86) V X (1,1) / Then we have x ' ( M ) in each of the domains (0,1) x (-L',0) and (0,1) x (0,L). The boundary conditions are Xtf ( l , l ) sin7rc/> (2.87) (2.88) (2.89) (2.90) xiu) r — * { 2 m ) in top and bottom domains correspondingly (with a2 = (x(l,l) + Ty ) /x ' ( l , 1) > 0). Using trigonometric formulas to simplify and substituting for <&j: [7(0, z) = (7(1,2) =0 Uz{<t>,-L) = Uz{<t>,L)=G. Seeking a separable solution in each of the domains we obtain by standard methods Ufa z)=(l- X H { 1 I \ ] ) ~ T Y ) (coshaz T tanhaLsinhaz) ^ §L(4>,z)=[l-XH(1,1) - Ty -1 + cosh a (X T z)\ sin7rc/> (2.92) X ' ( l , l ) J V coshaL y 7T Substituting this into the second continuity condition at the interface (2.85) and integrating: :(X(1,1) + TY) ( l - X h ( ^ ; y ) t a n h o L ] ' 9o + 9i = 7T" x ( i , i ) + ^y + ^ ] 1 [psin/?]2 • COS 7Tc/> • COS 7TC [St*x(l,l) + S t *Ty + pcos/?]2 Note, however, that from (2.76) the last term is O(e), and so as g\ is O(l) we have :(X(1,1) + TY) ( l - ^ ( ' j ^ ) t a n h a i ' (2.93) 9i -90 + X ( l , l ) + T y + ^ 20 12 1 - COS TT(f) (2.94) Chapter 2. Modelling of the flow using the relation (2.70) between x ' ( l , 1) and XH(1> 1) to simplify: 9i -90 + =F(X(1,1) + TY) ( l + ^ ^ ) t a n h a J ^ X ( l , l ) + ry + ^ l2 - COS 7T<j) Jl with the perturbed stream function in the moving frame of reference $i(<f>,z)= 1 + xm ( M ) Km - 1 + cosha(L=f z)\ sinrrci cosh OLL TT and as L —> oo-9i $ i ( < M -90 + 1 + T(X(1,1) + 7 V ) ( 1 - ^ ) ] " 7T X ( l , l ) + r y + ^ x m ( i , i ) 2 J 1 (—1 + cosh az ± sinh az - COS 7T(/) sinTf(/) (2.95) (2.96) (2.97) (2.98) In Figs. 2.2-2.4 we show main trends of these solutions, as we vary key physical parameters. In Fig. 2.2 we vary density of the displaced fluid pi and fix all other parameters at the values given. In Fig. 2.3 we vary the yield stress T y 2 and in Fig. 2.4 the consistency k2 of the displaced fluid. In each case the more elongated perturbation function g\ corresponds to the greater value of the varied parameter as we would expect. Note that for perturbation solution we require g\ = O(l ) << 1/e. A contour plot of the perturbed stream function $\(<j>,z) in the moving frame of reference is given in Fig. 2.5. 21 Chapter 2. Modelling of the flow -10.0 Figure 2.2: Plots of the interface perturbation function g\ for varying density of the displaced fluid p2- Least elongated to most elongated: p2 — 0.5, p2 = 0.6, p2 = 0.7, p2 = 0.8, p2 = 0.9, p2 = 1.0, p2 = 1.1. Other parameters fixed at: k\ = 1.0; fc2 = 0.9; m\ = m2 = 1.1; ryj = 1.0; 7V2 = 0.9; pi = 1; /? = 0; L = 10. Cn -10.0 Figure 2.3: Plots of the interface perturbation function g\ for varying yield stress of the displaced fluid Tyi- Least elongated to most elongated: T Y 2 = 0.6, Ty2 = 0.7, ry2 = 0.8, T Y 2 = 0.9, 7V2 = 1.0, TYI = 1-1- Other parameters fixed at: &i = 1.0; k2 = 0.9; m i = mo, = 1.1; r y i = 1.0; Pi = 1; P2 = 0.9; /? = 0; L = 10. 22 Chapter 2. Modelling of the how Cn - 1 0 . 0 Figure 2.4: Plots of the interface perturbation function g\ for varying consistency of the dis-placed fluid k2. Least elongated to most elongated: k2 = 0.6, k2 = 0.7, k2 = 0.8, k2 = 0.9, k2 = 1.0, k2 = 1.1. Other parameters fixed at: k\ = 0.9; m i = m2 — 1.1; TYI = 1.0; T y 2 = 0.9; Pi = 1; Pi = 0.9;-/? = 0; L = 10. 0 .4 0 .8 Figure 2.5: Contour plots of the perturbed stream function $i(c/>, z) in a moving frame of reference. Parameters fixed at: k\ — 0.8; k2 = 1; m i = m2 = 1.1; Ty i = 1.0; Ty2 = 0.9; pi = 1; p 2 = 0.9; 3 = 0;L=10. 23 Chapter 3 A n a l y t i c a l r e s u l t s In the previous chapter we derived the following differential equation for the stream function * : V a • S = - / , (3.1) where "~ T" n raTY S = faX( |V a t f |) raTY V a t f <==> \S\ > | V a * | H\Va*\ | V a * | = 0 <=>' | 5 | < H  : raTY H ' Subject to boundary conditions *(0 ,£ , i ) = 0 ,9* 0,0,*) = 0 tf(U,i) = Q(i), <9# and constitutive assumption 0 jjm+2 y K"l(m + 2) {x + ry/Hf ,m+l X + (m + 2)ry (m + l)H X < 0 , X > 0 , (3.2) (3.3) (3.4) (3.5) which implicitly defines x(|V a\IJ|). Here we show that there exists a solution to this equation and it is unique. For that we convert the problem to a minimization problem via a variational formulation, prove existence and uniqueness of a solution to the minimization problem and thus deduce the same for a weak solution of the original formulation. Thereafter we apply the augmented Lagrangian method to obtain an iterative algorithm that solves the variational formulation of the problem. 24 Chapter 3. Analytical results 3.1 Variational Formulation We note that the classical formulation (3.1 - 3.5) is not necessarily well-defined, due to the yield surface position x = 0 being initially undetermined. Therefore, we proceed to derive a formulation that is more open to analysis. The derivation is purely formal and we assume sufficient regularity of our solution and test function. So for the classical solution we assume * G CQ(O) and that * satisfies boundary conditions (3.4) with = (0,1) x (0,Z). First let us homogenize the boundary conditions by setting Note that can be determined by taking a linear combination of stream functions at £ = 0 and £ = Z. So in case of two fluids we can take = (1 — c)*o + c^z, where * a is the stream function at '£ = a and c is the concentration with c = 0 for the displacing fluid and c = 1 for the displaced fluid. Now let v, w G CQ(O) with w = v — u. Then from (3.1): # = u G CQ •02(fi). (3.6) wS • n ds (3.7) Substituting for w and S we have for |S | > Tyra/H: + ryrg ( V a * * + Vau) • {Vav - Vau) H | V a * * + Vau\ dCt + H (|va#* + vav\ - | v f l * * + vau|) dn (3.8) 25 Chapter 3. Analytical results Note that if \S\ < TyrA/H we have V a # = V a ^ * + Vau = 0, so S-Va(v-U) = S • ( V a r + V a « - (V a t t* + V 0 U ) ) = S • ( V f l * * + V a v ) < \S\\VA$>* + Vav\ < Z ^ | V a * * + V a w| . (3.9) Therefore (3.8) is valid for all S. Thus, a solution of (3.1 - 3.5) will also satisfy / | V a # * + V a u | - ( V a * * + V Q u) - ( V a « - V a u ) + ^ ( | V a * * + V a v | - | V a * * + V a u | ) + / ( 7 ; - w ) d Q > 0 (3.10) We now take (3.10) as the definition of our problem. In place of C2(Q) and C$(Q,) we assume that € V with M £ Vo (with VQ being the closed subspace of V with u = 0 on dVL). We determine the spaces V and Vb later. 3.J.l Immiscible fluids with no surface tension As shown in Chapter 2, if we wish to dispense with the concentration equation formulation, then the interface between two fluids is advected according to the averaged kinematic equation (2.55) and the following two conditions are satisfied at the interface between fluids 1 and 2 (respectively, the displacing and the displaced fluids): a* i a* 1 2 = 0 (3.11) l i.e. continuity of normal velocity. Also by continuity of pressure across the interface: [p]? = 0 (3-12) Additionally, ^ is assumed continuous across the interface. To derive the variational formulation in this case we proceed as above in section 2.1. So in each of the fluids instead of (3.1) we have: V a • S i = 0, for x e fii (3.13) V a • S2 = 0, for x e 0 2 , , (3.14) 26 Chapter 3. Analytical results with S i and S2 defined as in (3.2 - 3.5) in each of the fluids and Q,i,Q2 C ft are the regions where correspondingly fluid 1 or fluid 2 are present. Multiplying by a test function w € CQ(Q) and integrating over each subdomain: / S\ • Vaw dfti + / S2 • V a u ; dCl2 = <p wSi • n\ ds + f wS2-n2ds (3.15) corresponding to (3.7) in the variational formulation with the concentration equation above. Here the boundary integral terms on the right hand side vanish, except along the interface T, so we have f wSi-nids+f wS2-n2ds= / w(Si - S2) • n\ ds. (3.16) Substituting for S from derivation in Chapter 2: Sj = (raG^,-raG^) Thus where _ , dp p cos Brg dp rap sin B sin TT<J>\ ~ 1 r a ~ h i st* 'd<f> st* )j { 6 - U ) (S i - S 2 ) • m = ra [ti • V a ( p i - Pi)] - (fi - f2) • n i (3.18) * i = -ni,t) (3-19) is the tangent vector to the interface F and ^ = ^pco^ ^ ^ s M s i n ^ ^ ( 3 2 0 ) The first term in (3.18) is the difference in tangential derivatives of the pressure along the interface, which is zero by differentiating (3.12). Therefore, we have: J w(Si - S2) • ni ds = - J w ( / i - f2) • ni ds = — i> wfi • ni ds — f wf2 • n2 ds = ~ fi- V o ^ dfh - / f2 • Vaw d£l2 (3.21) 27 Chapter 3. Analytical results using integration by parts and since V a • fj = 0. Substituting this back into (3.16) and letting w = v — u: •r«x(|V atf* + V a u | ) , a . | V a # * + Vau\ - ( V a * * + Vau) • (VaV - V 0 U ) + ^ ( | V a * * + Vav\ - | V a # * + V a « | ) + /,• • Va(v - u) dQj > 0 H V u € C g , « e C g . (3.22) which is in fact equivalent to the variational formulation (3.10) since / = V 0 / and the term / f(v — u) dfi can be integrated by parts. Thus, working with either the concentration depen-dent formulation or considering the two fluid domains separately with continuity conditions on the interface as we have done here, leads to the same variational setting. The existence and uniqueness results together with the iterative solution algorithm that follow are thus applicable to both approaches. 3.1.2 Minimization Problem Consider the functional: v) = Jj I ^rjr- d s + ^ K * * + Va^ | - fv da (3.23) Equivalently, integrating the last term by parts and using the fact that v = 0 on SCI: J(v) = / y /  X\j2-L ^ + ^ | V a * * + Vav\ + f-Vav dfl, (3.24) where V a • / = / is as in (3.20). Let us show that the minimization problem J(u)< J(v),\/veV0, ueV0; (3.25) has a unique solution u G VQ, which solves the original problem in its variational formula-tion (3.10). 28 Chapter 3. Analytical results 3.1.3 E x i s t e n c e & Uniqueness Resu l t s Proposition 2.1. Let = / d s , r ( 3 2 6 ) with x as above. Then $ is strictly convex in R 2 . Proof. The Hessian matrix of is . &xx ®xv • H = ( " I , (3.27) with eigenvalues A i i 2 Ai,2 = \ + &vv) ± + * W ) 2 - 4($xx*yy ~ *ly)) ^ = ±(($xx + %y)±(($xx-%y)2 + 4$2xy)1/2) + X ' M r (4JZ4)2(^ - X'(r)) 2 + (^) 2(^ - X/(r)) x 2 + j / 2 r v / ; Kx2 + y2' K r = (^l+xir)^ ± *M_ x/( r) , where r = (x 2 + y 2 ) 1 / 2 . So A = 2 ^ ^ or 2%/(r) both of which are strictly positive. Thus $(x, y) is strictly convex in '. Proposition 2.2 J(v) as defined in (3.23) is strictly convex. Proof. We have for u,v G IR2; a € (0,1); 8 = 1 - a: -|Va**+Va(cm+/3u)|2 x ( s l / 2 ) ds dfi = ds dfi Jo + p [ r a [ ±1^-1 dsdQ, Jn Jo <a I ra I ds dfi In Jo s v . HVa**+V a^| 2 x ( s l / 2 ) 29 Chapter 3. Analytical results by Proposition 2.1. Also < f ^\Va** + Va{au + 0v)\dSl = Jo. n [ ^ | « ( V Q * * + Vau) + /?(V«tf* + Vav)\ dfi Jo. tt [ ^(a\(Va** + Vau)\ + 8\Va** + Vav\) dfi Jo. tt and, trivially, / /(cm + 8v) dD,< [ afu + Bfv dQ Jn Jo. Thus J(v) is a sum of two convex functions and a strictly convex function and is therefore itself strictly convex in R 2 . Now let us show that lim^v^+0OJ(v) = +00. Let K(x) = Jr-± j f d s + ^ \ x \ + f-(x-Vatf*) d n , (3.28) where x : K 2 -> R2 so that J(v) = K(VaV* + Vav). Proposition 2.3. K(x) —> 00 as Jn \x\1+1/m dQ, —> 00. Proof. From the constitutive relation (3.5) we have if x > TY/H: Hm+2 xm+l , (m + 2 )7y. V„\Tr = ± (Y+- i — M 1 a 1 km{m + 2)(x + rY/H)2^ {m + l)H! Hm+2 m+1 + 2 < (Y-\ '—Y) ~ km{m + 2) x2 m + lX' Xm(l + ^4) (3.29) fcm(m + 2) A ' v m + 1Also, for 0 < x < T y / ^ : | V * | < ^ ^ ( - V + . ( m + ^ -km(m + 2)2TYx/H H (m + l)H' Hm+2 m + 2 - 2fcm(m + 2 ) A v m + r < ^ Xmd \ m + 2) (3 30) ~ fc™(m + 2 ) X { L + m + lh [ 6 ' 6 [ ) ) 30 Chapter 3. Analytical results Thus, we have for % > 0 X > - a | V a ^ | 1 / r n , where Hm+2 , m + 2 A _ 1 / m n a = max | — ~ ^ ( H r) > 0. km(m + 2)K m + 1' Using this we have: Jn 1 Jo \X\2 „l/(2m) ,1/2 da - | / | | x | + / - V a * * dQ (3.31) (3.32) (3.33) where 5 = — f^f • Va<J/* dQ. Now, using Holder's inequality: - IVP f \f\\x\ dQ < c* / \x\p dQ Jn Un (3.34) I/P' where 1 < p < oo and c* = ^ / Q | / | p ' dQj with 1/p+l/p' = 1. Applying this with p = 1+1/m to (3.33): K{x) > & [ \x\1+1/m dQ - c* f \x\1+1^  Jn Un dQ l + l / n (3.35) where a = min( r a )a / ( l + 1/m) > 0. Rearranging, i f (a:) > a f Un \x\1+1/m dQ l + l / m /" | a ; | 1 + 1 / n Jn dQ + 5, (3.36) Thus K(x) —> oo as ( / n l a ; ^ 1 /™ d Q ) 1 / ( 1 + 1 / m ) -> oo. Theorem 2.4. J(v) —> oo as | | t » | | w i , i + i / m —> oo. Proof. First note that for x, y G R, with p = 1 + 1/m (1 < p < oo) and assuming 0 < x < y xp + yp< 2yp = 2(y2fl2 < 2(x2 + y2fl2 Thus, for a, b G L P (Q) we have H I L P + l i P < 2 | | ( a 2 + 6 2 ) 1 /2 | |P p 31 (3.37) (3.38) Chapter 3. Analytical results Now, by Poincare inequality we have that if w G L P (Q) and Vw G L p ( f i ) x LP{il) then | | W | | L P < i ^ | | V i y | | L P x L P , for some K > 0 (3.39) Thus, if v G W lv{Q) we have: IMIvpi.p = II^IILP + I I ^ ^ I I L P X L P < {l + K)\\Vv\\lPxLP < 2{l + K)\\{\Wv\)\\lP, (3.40) by above. Thus IMIiyi,i+i/m(n) -»• oo => | | ( | V v | ) | | i 1 + 1 / r o ( n) -» oo (3.41) Moreover, as ra is bounded and greater than zero we have IM I w i . i + i / m ( n ) ^ oo ^ | | ( | V a u | ) | | L l + 1 / m ( n ) oo (3.42) And as |V„** + Vav\ > \Vav\ - |V t ttf*|, (3.43) we have \ \ v \ \ w i , i + i / m { n ) -> oo =» HGVatf* + V a u|) |U 1 + 1 / m (n) -> oo. (3.44) And, finally by Proposition 2.3 we have \\v\\wi,i+i/m{n) * oo =y- tf(Vatf* + V 0 u ) = J(y) -+ oo, (3.45) as required. Summarizing, J(i>) is continuous, strictly convex and J(v) —• oo as I M I j y i . i + i / m ^ ) —> oo. Thus as W o M + 1 / m ( Q ) is a closed convex nonempty subset of W 1 , 1 + 1 / m ( S l ) , from standard opti-mization theory in Banach spaces it follows that the minimization problem (3.25) has a unique solution u G W 0 1 , 1 + 1 / m ( f i ) . Remark. In the preceding paragraph we in fact take m as the maximal m in Q. In case 32 Chapter 3. Analytical results of the concentration formulation this will be maxm(c(a:)), and for the immiscible fluids formulation it will be max{mi ,m2}. Moreover, writing J(v) = Jo(v) + Ji(v) where Mv) = J y J ^rj^1 ds - fv dfi (3.46) is strictly convex and Gateaux-differentiable on W 0 1 , 1 + 1 / m ( Q ) and Ji(v) = f ^ | V a * * + Va«| dfi (3.47) Jn tt is convex and continuous, we have (see Chapter V in Glowinski [17]) that the unique solution of the minimization problem is characterized by (Jo(u),v- u) + Ji(u) - Ji(V) > 0, Vt; G W o ' 1 + 1 / m ( f i ) , u G W o ' 1 + 1 / m ( 0 ) (3.48) where (JQ(V),W) is the Gateaux derivative of Jo: (Jo(v),w) = Jim-t->o t Jn | v a w + vav\ substituting into (3.48) we extract the variational formulation (3.10). Thus we conclude that the variational formulation has a unique solution in WQl'1+l/m(Q). 3.2 Solution by augmented Lagrangian method Above we have shown that the problem in its variational formulation has a unique solution. Here we consider a method for obtaining this solution through application of an iterative algo-rithm described in Glowinski [17]. For that we consider the minimization formulation of the variational problem defined above (see (3.25)): Minv€Vo{J{v)}, (3.50) 33 Chapter 3. Analytical results where J(v) = F(Vav) +G(v), with (3.51) F(q) = F0(q) + F^q), q £ H, where (3.52) Fo(q) = J-±J X^dsdn, qeH (3.53) Fi(q) = / ^ | V a # * + <?| dQ, q £ (3.54) G(v) = - f fv dn, v € V (3.55) Summarizing ,the results above we showed that Fo(q) is strictly convex and differentiable, Fi(q) is convex and continuous and G(v) is convex and continuous. Moreover, the min-imization problem has a unique solution in W 0 1 ' 1 + 1 / / m ( Q ) , so with VQ = WQ' 1 + 1 ^ m (Q) and H = L1+1/m(n) x L1+1'm{n). Note that for m < 1 we have that W 0 1 , 1 + 1 ^ m ( O ) is a subset of H^fl) which is a Hilbert space. However it is more common for the fluids used in the cementation process to be shear-thinning giving m = 1/n > 1. In that case we have WQ ,l+1^ m(n) 0 H^n) and so the solution doesn't necessarily lie in a Hilbert space. Application of the augmented Lagrangian method however requires the solution space to be a Hilbert space in order to obtain the convergence results. Nevertheless, once the variational problem has been approximated by some finite dimensional numerical method - finite element or, as we have done below, by a finite volume method, the approximate solution space Vh is finite dimensional - and thus is a Hilbert space. We hence take V = i?Q (Q) (so that Vh C V) and H = L2(n) x L 2 ( 0 ) , making it possible to apply the existence and convergence theorems to the iterative solution of the approximate problem. 3.2.1 Principle of the method We define a Lagrangian functional C associated with (3.51), by C(v,q,p)=F{q)+G(v) + {p,Vav-q), (3.56) 34 Chapter 3. Analytical results and for r > 0 an Augmented Lagrangian £ r defined by . £r(v,q,p) = £(v,q,p)+ ^\Vav - q\2 (3.57) Then, following Chapter VI in Glowinski [17] it is true that {u,p, A} is a saddle point of £ if and only if it is a saddle point of £r, Vr > 0. Moreover such u is a solution of (3.51) and p = S7au. Thus the original variational problem (3.10) is equivalent to finding saddle point of £ r which is accomplished by an application of an algorithm of Uzawa type. The resulting iterative algorithm is: \° e H given; then, A n known, we define un,pn,Xn+1 by £r(un, pn, Xn) < £(v, q, Xn), \/{v, q}eVxH, {un,pn} G V x H Xn+1 = \ n + pn{Vaun-pn), pn>0. The problem (3.59) is equivalent to solving two coupled variational inequalities: G(v) - G(un) + (A", Va(v - un)) + r (V . t t n - pn, V > - u11)) > 0, W G V > " G V , F(q) - F{pn) - (Xn, q-pn) + r(pn - Vaun, q~pn)> 0, Vq G H,pn G H The main drawback of this algorithm is that it requires solution of coupled variational inequal-ities at each step. To overcome this we uncouple the two inequalities in the natural way to obtain the following modified algorithm: p ^ A 1 G # given; (3.63) then p71"1, Xn known, we define un,pn, Xn+1 by 35 (3.58) (3.59) (3.60) (3.61) (3.62) Chapter 3. Analytical results G(v) - G{un) + (Xn, Va(v - un)) + r(Vaun - pn~\Va{v - un)) > 0, Mv £ V, un £ V, (3.64) F(q) - F(pn) -(Xn,q- pn) + r(pn - Vaun, q - pn) > 0, Vg € H, pn G H. (3.65) Xn+1 = x n + pn(yaUn _ pn^ ^ > g ( 3 6 6 ) In fact the first step (solving (3.64)) corresponds to minimizing Cr(v,pn~ l, Xn) with respect to v to get un and the second (solving (3.65)) to minimizing £r(un, q, A") with respect to q to obtain pn. 3.2.2 Application of the Algorithm In application to our particular problem with F and G as in (3.51 - 3.55) above we obtain: G(v)-G(un) + (Xn,Va(v-un)) + r(Vaun-pn-\Va(v-un)) • = [ -f(v - un) + Xn • Va(v - un) + r(Vaun - pn) • Va(v - un) dQ, (3.67) Integrating second and third terms by parts and noting that v, un = 0 on 8Cl we obtain G(v) - G(un) + (Xn, Va(v - un)) + r(Vaun - pn-\Va(v - un)) = / " ( - / - V a A n - r V a • (Vaun -pn))(v - un) dfi, (3.68) Thus the first step of the algorithm (3.64) is to solve f(f + V a A n + r V a • (Vaun - pn))(un - v) dfi > 0, v £ V, un £ V. (3.69) Jn This equation is satisfied by the solution to: rVa • Vaun = rVa • pn - VaXn - f, un £ V, (3.70) which is a Poisson equation on Cl, with the right hand side known at each iteration. 36 Chapter 3. Analytical results Going back to formulation of the second step of the algorithm we mentioned that finding pn at each step was in fact equivalent to minimizing Cr(un,q,Xn) with respect to q. So: pn = mi{Cr(un,q,Xn)} = m{{F(q) + G(un) + (Xn,Vaun-q)+r-\Vaun-q\2} +  rf\q\2-(Xn + rVaun)-q)c\n}, (3.71) Equation (3.71) is minimized when pn is given locally by: + y | V a # * + q\2 - (A™ + rVaun + r V Q * * ) • (q + V a * * ) } . (3.72) The expression above is minimized when (A n + rVaun + rVa**) is parallel to q + V a * * , since apart from the last term, the rest of the expression is a function of \q + V a * * | only, i.e. independent of its direction. So, letting q + VaV* = 9(Xn + rVaun + rVa**), (3.73) x = \Xn + rVaun + rV a.tf*| (3.74) we have to find the minimizer of M(6) = £ j f ( t e ) 2 ds + ?£\0\x + ^{Bx)2 - Ox (3.75) Now if x < s g * then ^\6\x-0x>0 (3.76) H and so with the rest of terms in (3.75) non-negative that implies 6 = 0 minimizes M, giving pn = - V a * * . If x > then we find 6 by setting ^ = 0, so raX(6x) + 2g£ + r ^ - x = 0 (3.77) which can be solved by. numerical inversion and thus giving p" = 8(Xn + rVaun + r V a * * ) - V Q * * . (3.78) 37 Chapter 3. Analytical results 3.2.3 Summary and convergence So, summarizing, the Augmented Lagrangian method applied to our particular variational problem results in the following algorithm: p°,\l eH given; (3.79) then p n _ 1 , A n known, we find un,pn, A n + 1 from rVa • Vaun = r V a • pn - V a A " - / , un G V, (3.80) f - V a * * • if x < 2t> p n = < (3.81) [ 0(A" + r V a u n + r V a * *) - V a # * if x > ^ A ^ 1 = A" + pn(S7aun - pn), P n > 0. (3.82) where x = |A" + r V a u " + r V Q * * | , and (3.83) TaXiOx) + ^ + rQx - x = 0. (3.84) Following Glowinski [17] this algorithm has the following convergence properties if 0 < pn < (1 + V5>/2: un -> u strongly in V (3.85) pn -* p strongly in H (3.86) An+i _ A n ^ o strongly in H . (3.87) where u is the solution of the minimization problem (3.23) and thus the variational formulation (3.10) for the approximate discretized problem. 38 Chapter 4 N u m e r i c a l R e s u l t s In the previous chapter we derived existence and uniqueness results for the model of the flow described in Chapter 1. Additionally, we applied the augmented Lagrangian method to the variational formulation of the equations for the stream function of the flow to obtain an itera-tive algorithm resulting in a sequence of functions converging to the solution of that formulation. In this chapter we outline the numerical solution of the equations of the algorithm and compare the results with the analytic solutions, in the case of concentric annulus, and with the perturba-tion solutions for small annular eccentricity. We apply a shock-capturing time-advance scheme to the concentration advance equation to arrive at a complete numerical simulation of the flow as described by the model given in Chapter 2. 4.1 Discretization and method of solution The algorithm resulting from the application of the augmented Lagrangian method to the weak formulation of the equations for the stream function of the flow is: P »°, A 1 G H given; (4.1) then pn 1,Xn known, we find un,pn, Xn+1 from ,n = rVa-pn-VaXn- /, un G V, (4.2) d(Xn + rVaun + rVa^*) - V a * * if x > 9g* - V a * * (4.3) 39 Chapter 4. Numerical Results = \ n + pn(Vaun-pn), P n > 0 . (4.4) where x = \Xn + rVaun + r V a * * | , and rax(Ox) + + r6x - x = 0, (4.5) (4.6) where V = H&(Q) and H = L2(fl) x L2(Cl), fl = (0,1) x (0,Z) and the relation between the modified pressure gradient % and the modulus of the stream function | * | is given from the constitutive assumption |V«* | = < 0 [ Km(m + 2) (x + ry/H)2 .m+1 x-(m + 2)TY X<0, X>0. (4.7) (m + l)H The result of solution of the algorithm is a sequence of functions un converging to u € HQ(Q) where the stream function * = \ F + u. Also we obtain pn £ L2(Q) x L2(Cl) converging to p = Vau from which the components of the velocity vector to be used in the concentration-advance equation can be easily calculated. We descretize the equations by 'applying a finite volume discretization with a staggered regular rectangular mesh, i.e we divide the domain Cl into Ni x Nj control volumes with each control volume being a rectangle of width A</> = 1/JVj in the ^-direction and height A£ = 1/Nj in the ^-direction. Denoting by CVij the control volume with its lower left-hand corner at coor-dinates ( iA0, jA£) we define the discretized functions v%, ^ at * n e corners of CVij and the discretized functions AJJ,p^ and concentration ch at its center - see Fig 4.1. This allows sim-plified calculation of gradients arising in the above equations as well as smaller spatial support of the numerical schemes for the chosen second degree of accuracy in each of the numerical calculations of the flux between the control volumes. We let the discretized concentration function ch be defined at the centre of each control volume and take values between 0 and 1 with ch = 0 in the displacing fluid (cement) and ch = 1 in the displaced fluid (mud). The density and rheological parameters are then also defined at the 40 Chapter 4. Numerical Results (iA4>! ((*+l/2)A0, (j +1/2) AO O jA£) O - A£, p£, c fe, V a ^ Figure 4.1: Schematic picture of the control volume CVij centers of the control volumes using linear mixing laws. So for example pfj = (l — cfj)po + cfjpi where po is the density of the cement and p\ is the density of mud. At the start of each run we first determine the "homogenizing" stream function Wh by set-ting = (1 — ch)^Q + Ch^z as in 2.1. Here * Q and ^ are the discretized stream functions at £ = 0 and £ = Z correspondingly. They in turn are determined from the boundary conditions which give the modified pressure gradient G = % + T y / H constant in the 0-direction at the top and bottom of the well. Thus, x at £ = 0 and £ = Z can be calculated by numerically inverting the' expression for the total rate of flow through the well. At each iteration of the algorithm we first solve the u 1^  advance equation (4.2). We applied the standard point Gauss-Seidel method with second-order accurate spatial fluxes to solve this modified Poisson's equation. This seems especially appropriate as v% converges after some it-erations of the algorithm and so taking as starting point of Gauss-Seidel iterations each time improves the convergence time. The same is true after application of the time-advance scheme: we then take the final value of v% at the previous time-step as the starting point for Gauss-Seidel iterations. Over-relaxation was also used to improve the speed of convergence even further in individual cases. 41 Chapter 4. Numerical Results The second step in the algorithm is the advance of p\. We achieve this by first calculating x from (4.5), determining if the fluid is yielded at that point and if so performing the numerical inversion in (4.6) to determine 0 and thus p\. Note that (4.6) can be differentiated analytically and a Newton-Raphson method used in its inversion converges in few iterations to the desired tolerance. 4.2 Convergence We take the norm of Ap\ =pl+l-pl as the measure of convergence of the numerical algorithm. This appears appropriate as it is the velocity field obtained directly from p\ that is used in the subsequent concentration advance step. Here we consider two norms of Ap 1^ : L 2 (Q) and L°°(0). While the first provides an adequate measure of convergence to the solution over the whole domain, the second gives a better impression about the convergence near the interface and yield boundaries - where it is the slowest. The plots given in Fig.(4.2) are of the logarithm base 10 of the I? and L°° norms of Ap% versus the number of iterations together with surface plots of a component of Ap^ after 100 iterations. Different interface configurations and situations are presented. It can be noted that absence of a yield boundary within the flow (so that both of the fluids are fully yielded in the annulus) and lower eccentricity improve the convergence dramatically. In Fig.(4.2) this can be seen from the difference in the convergence histories between cases a) and b). The rheological parameters, densities and the interfaces between fluids are the same in both cases - however b) has a higher eccentricity. From the AJJ-advance equation (4.4) we see that as AJJ converges, we should have p\ —> V a w ^ . The convergence history of the L 2 and L°° norms of APU = p\ — V a w ^ = ( A ^ + 1 — A^) /p n for the same set of conditions as for A p \ above, is shown in Fig.(4.3). The convergence mirrors that of ApJJ with the notable exception of the behavior in the unyielded region (this can be seen in case c). The reason for this seems to be that whereas in this algorithm vfi converges to 42 Chapter 4. Numerical Results Figure 4.2: Plot of Logio of the l? (solid line) and L°° (broken line) norms of ApJJ against the number of iterations together with surface plots of the ^-component of Ap^ after 100 iterations (top to bottom): a) horizontal interface with e = 0.1, b) horizontal interface with e = 0.6, c) slanted interface with e = 0.8 and high yield stress of the displaced fluid. The mesh size here as well as in all of the subsequent examples was taken to be 20 x 80. 43 Chapter 4. Numerical Results > a modified stream function and p\ converges to the modified velocity field correspondingly, Xfc converges to the modified pressure gradient field by virtue of the u^-advance equation (4.2) and in our model the pressure gradient is undefined in the unyielded region. This implies that in the "unyielded" region (i.e. where |A£ + rVau^ + r V a \ ] / * J < Tyra/H) at each step of the algo-rithm we set p\ = - V B * j , solve V a • V a W / J = V a - p \ = - V a • V a * ^ and advance as before. Lastly we show that the original field equations are satisfied by the iterative solution with an increasing degree of accuracy. For that we set = p1^ + *£ and substitute into the field equations (3.1 - 3.3) to obtain S^. Below in Fig. (4.4) are the plots of the L2 and L°° norms of the residual Res = V a • + fh with fh being the discretisation of the function appearing on the right-hand side of the field equation as defined in (2.48). The geometric, rheological and density parameters are the same in each case as in the plots for APU and Apfc. The residual .decreases in the fashion similar to the convergence of the algorithm's functions examined above. It should be noted that in the more extreme case c) where the eccentricity and yield stresses are set high and the yield boundary runs through the domain, the residual plotted is large in the control volumes adjacent to the yield boundary. This is a numerical artifact arising from the fact that in order to compute V a • numerically second order accurate in space from the functions in the algorithm, we take a larger spatial support than we do to calculate fh-Thus, near the yield boundary part of the support for V a • S1^ lies inside the unyielded part of the domain where it does not converge to — fh by the virtue of the underlying model. This should not be a problem for the overall solution as we feed back only the velocity field into the concentration advance equation which, in turn, is extracted from p\ that has the same spatial support as fh and converges numerically on the entire domain as we illustrated above in Fig. (4.2). The norms of the residual Res in this case are calculated over the control volumes away from the yield boundary. 44 Chapter 4. Numerical Results N u m b e r o f I t e r a t i o n s V 5 4 Figure 4.3: Plot of Logic, of the L2 (solid line) and L°° (broken line) norms of APU against the number of iterations together with surface plots of the (/(-component of APU after 100 iterations (top to bottom): a) horizontal interface with e = 0.1, b) horizontal interface with e = 0.6, c) slanted interface with e = 0.8 and high yield stress of the displaced fluid. 45 Chapter 4. Numerical Results 14 Figure 4.4: Plot of Logw of the L2 (solid line) and L°° (broken line) norms of the filed equations residual Res = V a • + fh against the number of iterations together with surface plots after 100 iterations (top to bottom): a) horizontal interface with e = 0.1, b) horizontal interface with e = 0.6, c) slanted interface with e = 0.8 and high yield stress of the displaced fluid. 46 Chapter 4. Numerical Results 4.3 Time advance scheme The concentration advection equation | [Hrack] + ^ [Hv cfc] + ^ [Hraw ck] = 0, (4.8) is used to calculate the advance of the concentration and thus the interface after solving the field equations and extracting the velocity field (v, w). To simplify the notation we rewrite this equation as: 8U dF OG n , . -bJ + ^ + ^ = ° ( 4 - 9 ) where U = Hrack F = Hvck G — Hraw ck (4.10) This is a hyperbolic conservation law and different methods for obtaining its numerical solution exist. In our specific situation we have a large concentration gradient at the interface between the fluids so a time advance scheme must be chosen which captures propagation of this interface accurately. High order time advance schemes (second order accurate and above) such as Lax Wendroff are known to suffer numerical dispersion near discontinuities and regions of high gradient. This results in unphysical oscillations producing under- and overshoots in the conserved quantity U near the interface. Conversely low order order schemes (first order accurate) such as donor cell or upwind Euler don't produce dispersion but, because of their low accuracy, numerical diffu-sion is introduced. This has an effect of "smoothing out" the interface which is much greater than is present physically. The schemes which effectively minimize both the numerical diffusion and dispersion by hy-bridizing upwind low-order schemes and higher order schemes while limiting the total numerical flux functions F and G are known as Total Variation Diminishing (TVD) and Flux Corrected Transport (FCT) schemes. They have been used extensively for the solution of the hyperbolic conservation problems such as those of shock and signal propagation. 47 Chapter 4. Numerical Results For the examples following we used an F C T scheme with donor-cell upwind discretisation for the low order scheme and central difference discretisation for high order scheme. This was found to produce good results in previous applications to the concentration advance equation above and is used currently in industrial applications. A two-step MacCormack T V D scheme with different limiters was also tried out and found to produce similar results varying slightly with the types of limiters used. A brief outline of the F C T scheme used is as follows. First, low order fluxes FL and GL and high order fluxes FH and GH are computed using a donor-cell upwind and central differ-ence discretisations correspondingly. From those we compute antidiffusive fluxes A1 = FH — FL, Ai = GH — GL and the low-order time advance solution Utd from the low-order fluxes FL and GL. Secondly, we limit the antidiffusive fluxes by setting AiC = Ai C* and AiC = A1 CK Here the limiting functions 0 < C\ C J < 1 are chosen to limit U(t + At) so that it does not exceed some maximum Umax(t) or fall below some minimum Umm(t). In turn we choose Umax(t) to inhibit the overshoots in ck as the maximum value of U(t) over the neighboring cells and analogously for Umm(t). Lastly, we apply the limited antidiffusive fluxes to the low order time advanced solution to get the final time advanced solution: U(t + At) = U t d - - ^ \AiC - A ^ c + Ai° - A ^ " 1 * 0 ] (4.11) where AV = A0A£. The procedure used for each of the following simulations was to first set up an initial con-centration field (usually corresponding to a horizontal interface in the middle of domain) and then run the augmented Lagrangian solver followed by the time advance scheme above a desired number of times. To limit the size of domain (and thus drastically reduce the computational time) we "followed the interface" by shifting the concentration field down each time after a certain number of time steps by effectively adding a fresh row of control cells at the top of the domain and removing a row from the bottom. The number of time steps between each 48 Chapter 4. Numerical Results adjustment was chosen so as match the average speed of the flow. This has the effect of always placing the interface near the middle of the domain and allowing us to restrict the vertical size of computational domain to just a few units either side of the interface. 4.4 Simulation results Here we give several results from the simulation of the displacement flow using the ideas de-scribed above for solving the flow equation and the concentration advection equation. For the record, we use a 20 x 80 mesh with rectangular control volumes of width A(f> = 0.05 in the (^-direction and of height A£ = 0.1 in the ^-direction. We pick the length of a time-step A i from the Courant-Friedrichs-Lewy (CFL) stability condition: A t = — ~~i > (4-12) -^-V -v- -1-W vmax T w max where Co(< 1) is the Courant number and Vmax, Wmax are the maximum absolute values of numerical velocities in the <f> and £ directions respectively extracted from the the solution of the <£-field equation. A common procedure in such cases is fixing Co and adjusting the time step after each time advance accordingly. However it can be noted that in our case both Vmax and Wmax vary little (typically by less than 5%) throughout the length of the simulation. Thus we fix A t based on the above formula at the beginning of each simulation run which addition-ally allows us to easily " follow the interface" by periodically shifting the concentration field as described above. Typical values for a medium-eccentricity "unstable" displacement situation were Vmax = 0.5, Wmax = 4 so with Co = 0.5 we would take A i = 0.01 in this case. First we examine the evolution of the interface in the case of a concentric annulus with e = 0. The analytical solution for the interface in this case was'obtained in 2.3.2. Here, as in most of the simulations following, we set the interface to be initially horizontal and positioned in the middle of computational domain. The results for three different situations are given in Fig. 4.5. Plot of the projected interface position after 30 time units is given together with the time track of the interface position on the side 0 = 0. The "interface" is taken to be the contour .c^ = 0.5. 49 Chapter 4. Numerical Results The situations considered are those described in Chapter 1. The first is that of a vertical well where the interface stays horizontal. The second is for a horizontal well, where the steady-state interface is given by g(<f>) = - ^(/j+Vy]2 cos(7T0)/(7r5t*). The third is for an inclined well (8 = 5TT/12 was taken) with fluids of identical rheologies but a non-zero density difference, where the steady-state interface is given by g(<f>) = - tan/? cos (TT<P)/T. In all cases the simu-lated interface converges to that predicted analytically, they coincide on the plots shown and the average of the position of the interface on the wide side at the end of each simulation agrees with analytical result. Note that from the mass conservation considerations if the interface on the wide side becomes steady then the interface on the narrow side is also steady. Secondly, we compare the results of the flow simulation in case of small annular eccentricity to those obtained though perturbation methods in Chapter 1 with the formula for perturbed interface given by (2.95). Numerical simulation of such small perturbations is somewhat prob-lematic as the smaller mesh size we use to achieve sufficient resolution in the ^-direction not only increases the overall grid size but also necessitates a smaller time step by the C F L con-dition (4.12). Results for two different sets of parameters are given in Fig. 4.6 together with track of the simulated interface position (in frame of reference moving with the average speed of the flow) on the wide side. In the flow simulations we start with a horizontal interface as before and the mesh size is 80 x 20 with grid spacing of 0.05 in both 4> and £ directions. The computed solution at t = 30 is within 0.02 of the analytic solution in each case, note that the vertical scale on Fig. 4.6 is smaller than that in Fig. 4.5. 50 Chapter 4. Numerical Results 6 . 0 5.5' 5 . 0 ' J . 5 o II 4 . 0 ' '5 3 . 5 ' 3 . 0 2 . 5 0.6 1 . 0 25 30 10 15 0.2 0 . 4 0 . 6 0.8 1. 0 10 15 20 25 30 35 40 Figure 4.5: Plot of interface in a concentric annulus at t = 30 together with interface track at 0 = 0 (top to bottom): a) vertical well /3 = 0 (e = 0; p\ = 1; p 2 = 0.6; T y i = 1; T y 2 = 0.8; k\ = 1; A;2 = 0.8), b) horizontal well /? = n/2 (e = 0; p\ = 1; p 2 = 0.6; T y i = 1; T y 2 = 0.8; fci = 1; fc2 = 0.8), c) inclined well with angle of inclination 8 = 5?T/12 and identical rheologies (e = 0; pi = 1; p 2 = 0.8; T y i = 1; T y 2 = 1; k\ = 1; fc2 = 1)- Inverse power index mi = m 2 = 1.0 throughout. 51 Chapter 4. Numerical Results . 0 0 • 3 . 0 Figure 4.6: Plot of the interface position given by the flow simulation (solid line) and pertur-bation solution (broken line) together with the track of the interface position on the wide side of the annulus (top to bottom): a) e = 0.1; p\ = 1; pi — 0.8; ry\ = 1; ry2 = 0.9; k\ = 1; k2 = 0.9, b) e = 0.2; pi = 1; p 2 = 0.8; Tyi = 1; ry2 = 0.8; k\ = 1; fc2 = 0.8. Inverse power index m i = = 1.1 and angle on inclination /3 = 0 throughout. 52 Chapter 5 L u b r i c a t i o n m o d e l One of the main objectives of modelling the annular flow considered here is to simulate behavior of the interface between fluids so as to determine the type of displacement that occurs given the specific fluid rheologies, densities and geometry of the well. From physical considerations and numerical simulations of the flow, such as those described above and in [15] it is evident that one of three main situations can occur: 1. a " steady displacement" with both fluids fully yielded and the interface stationary in the frame of reference moving with the average speed of the flow; 2. an "unsteady displacement" with both fluids fully yielded and the displacing fluid shooting up on the wide side of the annulus with the interface elongating; 3. a "partially unyielded displacement" with either the displaced fluid or both fluids un-yielded (stationary) on the narrow side of the annulus and the interface between them elongating. In previous chapters we developed an iterative numerical method of solution for the flow model and showed its functionality. Although the computational speed of the resulting simulations is much faster than that of a fully 3-D flow simulation, it should be noted that for larger domains with smaller mesh spacings, highly eccentric annuli or long simulation times, the computational time can be significant. This could be improved for example by employing faster algorithms for solution of the Poisson equation for advance of un above or a more extensive exploitation of the over-relaxation constants pn and r. Nevertheless, solution of the 2-D model is slow and it would be advantageous to determine which of the three types of displacement occurs in a 53 Chapter 5. Lubrication model particular situation without the lengthly process of simulating the 2-D flow. In this chapter we employ techniques used in modelling of lubrication and thin film flows to derive a criteria for determination of the type of annular displacement occurring and compare these with results of numerical simulations performed using the algorithm for the. 2-D flow, described in the previous chapter. 5.1 M o d e l D e r i v a t i o n As a starting point we take the model equations for immiscible fluids with no surface tension described in Chapter 2. This is given by the field equations (2.52 - 2.53) V a • Si = 0, for x e flu (5.1) continuity conditions (2.54, 2.56 & 2.57) on the interface £ = h((f>,t): bl? = 0; = 0, (5.2) which lead to: [* • Vop]?"= [t • V a * ] ? l_dp dpdh12 ra dcf) <9£ d4> J _ d * dVdh ra d4> d£ dcj) = 0 (5.3) I 2 = 0 (5.4) with t = (1, | | ) being tangent to the interface and the kinematic condition at the interface (2.55): dh v dh dt rad(j) The overall idea here is to assume a highly elongated interface which results in several simplifi-cations to these equations. From there we can deduce that if, based on the simplified equations, the interface elongates further it indicates that the situation would lead to an unsteady dis-placement to start off with. Whereas if the interface "levels off" it would indicate that a steady displacement would normally take place. Finally, we are also able to determine whether or not a mud layer (fluid 2) on the narrow side would be mobile in the simplified model. Thus, we can distinguish the 3 cases listed in the beginning of this chapter. 54 Chapter 5. Lubrication model 5.1.1 Lubrication assumptions In order to apply the lubrication approximation we assume the following: Streamlines are pseudo-parallel so that the main velocity component is in the ^-direction. Thus we have «W (5"6) Interface is pseudo-parallel i.e. is highly elongated in the ^ -direction. Denoting the inter-face by (j) = this translates into | | | < 1 . (5-7) Modified pressure gradient S has its main component in' the ^-direction: \S(-\< 1^ 1, (5.8) as follows directly from the first assumption and the definition of S in the case the fluids are yielded. When either of fluids is unyielded this is in fact an additional assumption. Intermediate scaling. In the derivation of the flow model we assumed that although the angle of inclination of the well 8 and the mean half-gap ra vary with £ though the length of the well, they do so gradually and on a length scale of many annular diameters. When performing most numerical simulations here we also assume for simplicity that 8 and ra are fixed in the computational domain. This would correspond to considering a length of the well a few annular diameters long where they can be assumed to be constant. Here we use a similar argument to fix 8 and ra. To make it more rigorous, introduce an "intermediate" length scaling factor e such that rrr* and e << 1 where in the notation of Chapter 2, £tz — ^ is the total dimensional length of the well and r* is the global average mean radius of the annulus. Typical values for these are 55 Chapter 5. Lubrication model itz — £bh ~ 500 meters and rrr* « 0.3 meters so we consider e « 10 to be a suitable intermediate length scale allowing that. 5.1.2 Scalings and approximation derivation Using the intermediate scale described above we re-scale the axial and time variables by setting ( = eZ; t = et (5.10) also set W = w; V = ev; P = ep (5.11) for the components of the velocity vector and the pressure correspondingly and write *(</>, £) for the stream function as before. For C = 0(1) we have by the argument above, 3 is constant, ra = 1 and H = l + ecos7rc/>. Under these scalings the kinematic equation (5.5) for the interface 4> = <&(C>*) becomes: % + W%-V ( M 2 ) at the interface. Using the lubrication assumptions for the stream function * and modified pressure gradient S, the field equations (5.1) become: &(s*O)=0 ("3) to leading order in each fluid. Thus is constant in each fluid domain and so from the derivation of the field equations given in Chapter 2 we have _ dP pcosd * = ~~cX~^F= • ( } and by lubrication assumption 1 dP psm6smir<t> S< = -edj + St* < < X ( 5 ' 1 5 ) thus ^ = ep sin 0 sin IT(f)/St* to leading order. Now consider the continuity condition (5.3) for the tangential derivative of p at the interface: dP_ dP d(f>il2 = 0. (5.16) l • 56 Chapter 5. Lubrication model Thus, as 4? = 0(e) from above we have dP = 0 (5.17) at the interface. Thus ^ = Pc(C) is constant for fixed £ for all </> 6 [0,1]. From the lubrication assumption we have • \G\ = \S\ = \S+\ = \P< + p cos (3. St* ' (5.18) to highest order in each fluid. We can find P^(C) from inverting the flow-rate equation I ^< G '» d * + lw ( G 2 ) d * = 1 (5-19) for fixed £ and a unit flow rate through the annulus, where Pi cos P. G i = |P< + G2 = \Pr + St*  1 P2 cos (3 St*  1 (5.20) (5.21) and the relation between G and | | - is given from the constitutive assumptions by function F(G) as in (2.40). Thereafter we can obtain the stream function * from P^ by integrating ^ = F(G) with respect to (f> and using the boundary condition *(0 = 0) = 0. From there we can obtain the velocity field W _ T 0 * Tllty' V HOC (5.22) 5.2 Kinematic equation and interface stability Now we examine the evolution of the elongated interface from the approximation derived. The kinematic condition at the interface (5.12) gives Rewriting first term and grouping the second and third together we have: (5.23) d_ dt r4>i / H(4>) Jo dcf) (5.24) 57 Chapter 5. Lubrication model write $i(C,t) = / H(4>)dct>, (5.25) Jo g($4) = ViMMCt))) (5.26) Using this notation we obtain a hyperbolic equation for propagation of the "interface": f + ! , ( * , ) = », (5.27) with boundary conditions q(0) =0; 9(1) = 1- (5-28) Nota that <I>; is actually a "volumetric" position of the interface, i.e. the volume fraction of fluid 1 at depth £. From this we obtain the speed of propagation of the interface in the (- direction: d^dq dq_ 1 d^dC d * i ' K ' provided no shocks occur. Due to hyperbolicity of (5.27) shocks can occur. They correspond to patches of the interface where | ^ | becomes large and the interface becomes near-azimuthal. In order to calculate the speed of propagation of the interface at such a shock we use mass conservation considerations (effectively Rankine-Huguenot conditions) . Dividing shocks into three types depending on their azimuthal position in the domain: Shock on wide side occurring for $ G [0, $i,w]. Denoting the speed of propagation of the interface in the (-direction as WitW we should have the following holding true for the flux function q: q(* = **,«) = $i,wWi,w (5.30) • so as to the right of the shock we have Wj($ = $i ) l u) = q'{$i,w) and Wj should be continuous we have: WitW = q'($i,w), where (5.31) q(*i,w) = $i,wq'($i,w)- (5.32) 58 Chapter 5. Lubrication model Shock i n central region occurring for $ G [$+"c, ^  ]. Denoting the speed of propagation in the ^-direction by Wj ) C we have analogously to the above: g(* = <E>+ ) - <?(* = ) = ($+ - * r (5.33) and also wi,c = 9 , « c ) = 9'(<&r;c). (5.34) Shock on narrow side occurring for $ G 1]. Denoting the speed of propagation by Witn we obtain as before 5.2.1 Determining the displacement type Now we are able to classify the displacements according to the behavior of the flux function - which is the stream function at the interface. Our classification will be based on the wide and narrow side interface velocities, which me must determine from q. The behavior is complicated by the presence of shocks. After computing q for given fluid densities, rheologies, angle of inclination and eccentricity, we can conclude the following: If q'(l) = 0 then a mud channel occurs on the narrow side of the annulus and the displaced fluid is unyielded there. The displacement is thus not steady and the interface continues to elongate. A graph of such q with its derivative is given in Fig.5.1 (a). If q'(l) ^ 0 then the fluid is yielded everywhere. In this case we first have to check if any shocks such as described above occur and determine their location. No shocks occur if q" is of one sign - so it is either strictly positive or strictly negative ev-erywhere. If q" > 0 everywhere (Fig.5.1 (b)), then q' is increasing and q is strictly convex then the speed of propagation of the interface in the (^-direction on the wide side of the annulus is WitW = q'{0) and is less then the speed at the the narrow side W^n = c/(l). Thus the interface Wi,n = q'($i,n), where (5.35) 1 - ? ($ i ,n ) = </($i,n)(l - $ i ,n ) . (5.36) 59 Chapter 5. Lubrication model becomes less elongated in the ^-direction and we can conclude that the displacement would be steady. If q" < 0 everywhere (Fig.5.1 (c)), then q' is decreasing and q is strictly concave then WiiTl = q'(0) > q'(l) = Witn- Thus in this case the interface comes increasingly elongated and we can conclude that the displacement would be unsteady with displaced fluid slumping on the narrow side and the displacing fluid "shooting up" on the wide side. Shocks do occur if q" — 0 for some $ € (0,1) (Fig.5.2 (a)&(b)). We determine their po-sition and boundaries from (5.32) & (5.36) and speed of propagation from (5.31) & (5.35). Comparing speed of the shocks on the wide and narrow side, if present, we conclude that the displacement is steady if W^w < W^n and unsteady otherwise. One other situation of note that can arise is that of a slumping flow - if the density differ-ence between the two fluids is too great then we can have the stream function q decreasing for some values of $ i (Fig 5.3). This corresponds to parts of the interface moving in the negative ^-direction opposite the average velocity of the flow due to significant buoyancy effects. For example in case of the displacing fluid having a much greater density than the displaced fluid we observe q'($i) < 0 for small - thus the interface moves downwards on the wide side of the annulus, indicating a stable displacement as we would expect. To summarize, the behavior of q for different rheologies and physical parameters is quite rich. We have shown only a selection of the potential behavior. 5.2.2 Numerical results Here we provide some results numerically obtained from the lubrication model. Apart from de-termining the type of fluid displacement in a particular situation with given rheologies, densities and well configuration, it is of interest to examine the change in the displacement type as we vary several parameters. To demonstrate the trends we plot the differential velocity Wj i W — W j i n at the sides of the annulus together with the displacement type we thus expect for a range of varying parameters. 60 Chapter 5. Lubrication model Figure 5.1: Plots of the stream function q and its derivative with respect to <& versus the interface position 4> for different situations. Top to bottom: a) mud channel on the narrow side (e = 0.8; p\ = p2 = 1; TYI — 1; T y 2 — 1-8; k\ = 0.8; fc2 = 1); b) steady displacement, no shocks (e = 0.1; p\ — 1; p2 = 0.8; T y i = 1; TYI = 0.8; k\ = 1; fc2 = 0.8); c) unsteady displacement, no shocks (e = 0.1; p\ = 0.8; p 2 = 1; T y i = 0.8; T y 2 = 1; k\ = 0.8; fc2 = 1). Power index mi = m 2 =4.0, angle of inclination [3 = 0 throughout. 61 Chapter 5. Lubrication model Figure 5.2: Plots of the stream function q and its derivative with respect to <& versus the interface position <f> for different situations. Top to bottom: a) shock on the wide side (e = 0.6; px = l ; p2 = 0.8; T y i = 1; Ty2 •= 0.8; fo = 1; = 0.8); b) shocks on wide and narrow sides (e = 0.1; P l = 0.9; p2 = 1; T y i = 1.2; T y 2 = 0.8; fo = 1.2; fo = 0.8). Power index m i = m2 = 1.0, angle of inclination /3 = 0 throughout. - i . o Figure 5.3: Plots of the stream function q and its derivative with respect to $ versus the interface position 4> for a slumping flow situation (e = 0.6; p\ = 1.2; p2 = 0.6; T y i = 1.2; T Y 2 = 0.6; fo = 1.2; fo = 0.6). Power index m i = m2 = 1.0, angle of inclination /3 = 0. 62 Chapter 5. Lubrication model 0.7 0 . 8 1.0 1.1 1.2 2 2 0 Figure 5.4: Contour and relief plots of the differential velocity Wi<w — Wi,n for varying density p2 and yield stress Ty2 of the displaced fluid and other parameters fixed (given in text above). On the left plot, darker shaded region corresponds to predicted unsteady flow and lighter to predicted steady flow. For Fig. 5.4 we fix the angle of inclination j3 = 0, the eccentricity e = 0.2, the inverse power index of both fluids mi = m2 = 1.1, their consistencies fo = 1.1, fo = 0.8, density of the displacing fluid p\ = 1 and yield stress of the displacing fluid Tyi = 1.1. We vary the density of the displaced fluid p2 between 0.6 and 1.2 and yield stress of the displaced fluid ry2 between 0.6 and 1.4. As expected we see that with increasing p 2 and ry2 the flow becomes unsteady. For Fig. 5.5 we now vary the density of the displaced fluid p2 between 0.5 and 0.9 and consis-tency of the displaced fluid fo between 0.4 and 1.2. We fix the angle of inclination 8 = 0, the eccentricity e = 0.3, the inverse power index of both fluids mi = m 2 = 1.1, their yield stresses r y i = 1.2, r y 2 = 0.8, density of the displacing fluid p\ = 1 and consistency of the displacing fluid fo = 1.2. Here the flow is slumping for smaller values of p\, fo, steady for intermediate values and unsteady for larger values as expected. Lastly, illustrated in Fig. 5.6 we vary eccentricity e between 0.2 and 0.8 and density of the displaced fluid p2 between 0.4 and 1. The other parameters are fixed at following values: fo = 1.0; fo = 0.8; m i = m.2 = 1.1; Tyi = 1-0; r y 2 = 1.2; pi = 1; 8 = 0. The predicted flow is steady for smaller eccentricity and density of displaced fluid, unsteady for intermediate values 63 Chapter 5. Lubrication model k2 p2 Figure 5.5: Contour and relief plots of the differential velocity Wi ) 1 f l — Wi, n for varying density P2 and consistency k2 of the displaced fluid and other parameters fixed (given in text above). On the left plot the darkest shade corresponds to predicted unsteady displacement, the lighter shade to predicted steady displacement and the lightest shade to predicted slumping (steady) displacement. with a mud channel on the narrow side predicted for larger values. 5.2.3 C o m p a r i s o n w i t h flow s imula t ions Now we examine how the results obtained from the lubrication approximation compare with those obtained by simulating the full 2-D flow using the methods described in previous chapters. First we examine how the prediction of a steady or unsteady flow compares with the results produced by flow simulation. Here, we fix all parameters but one and examine where transition from steady to unsteady displacement occurs as predicted by each model with variation of that parameter. In Fig. 5.7 we vary p2 and fix all other parameters at values specified. Plot of differential velocity WiiW — Wi^n calculated from the lubrication model is given together with results of 2-D flow simulation at different values of p2. In the lubrication model we consider negative value of WitW - W\,n to indicate a steady displacement and positive an unsteady displacement. Each of the 2-D simulations was done on a 16 x 100 mesh with grid spacing 0.0625 and 0.1 in the 4> and £ direction, respectively, with initial interface azimuthal. Determination of whether 64 Chapter 5. Lubrication model Figure 5.6: Contour and relief plots of the differential velocity W^w — Wi, n for varying density of displaced fluid p2 and eccentricity e of the annulus and other parameters fixed (given in text above). On the left plot the black region corresponds to predicted mud channel on the narrow, the lighter shade to predicted unsteady displacement with fluid yielded on the narrow side and the lightest shade to predicted steady displacement. displacement is steady or unsteady in a particular simulation run was done by examining the interface track on the wide side in a frame of reference moving with the average speed of the flow. There are however problems with such determination as shown below taking the particu-lar set of parameters given as an example. For small values of p2 - (0.4 or 0.5) the interface track revealed that the interface became fully developed within few time units with the position of the interface at the wide side in the frame moving with average speed of the fluid g(<f> = 0,t) asymptotically approaching a value few units above the starting level. As we increase p2 we can see that the "development time" becomes larger with increasing p2 and limiting value which g(0,t approaches becomes larger. In order to conclude whether the displacement is steady or not from the flow simulation we have to run the simulation for at least the length of this development time. However as this time gets increasingly greater and we have to take a larger computational domain in the £ direction to accommodate for the greater value to which g{0,t) approaches, the computational time increases greatly for larger values of p2. This is best illustrated in Fig. 5.9 d) where we set p2 = 0.7 and track the position of the interface for 40 time units (4000 time steps) with g(0, t) 65 Chapter 5. Lubrication model c • H s I s - H s 2 . 0 1. 5 1. 0 0 . 5 - 0 . 5 - 1 . 0 -1 . 5 - 2 . 0 5.00 4.75 4.50~ "3.75 3.50" 3.25" ^^^^ 5.00 4.75 4.50 .25 .00' 0-3.75 3.50" 3.25 Figure 5.7: Top: plot of the differential velocity WiiW — Wi>n calculated using the lubrication model versus the density of the displaced fluid pi with other parameters fixed (e = 0.4; k\ = 1.0; ki — 0.8; m i = mi = 1.0; TYI = 1-0; Tyi = 0.8; p\ = 1; /3 = 0). On the plot, solid circles indicate the values of density for which 2-D flow simulation indicated an unsteady displacement and empty circles where a steady displacement was achieved. Below are plots of the interface position at the wide side of the annulus (4> = 0) from the simulated flow in frame of reference moving with the average speed flow for (to to bottom, left to right): a) pi — 0.4; b) pi = 0.5; c) pi = 0.6; d) pi = 0.7; e) p2 = 1.0; f) p2 = 1.1 66 Chapter 5. Lubrication model Figure 5.8: Right: plot of the interface position at the wide side of the annulus from the simulated flow in frame of reference moving with the average speed of the flow for p2 = 0.7 with interface evolving from an initially elongated configuration (plotted on the left). Other parameters fixed at the same values as in Fig. 5.7. steadily increasing at a nearly constant rate as shown. We cannot however conclude that the displacement is thus unsteady as running the simulation for 80 time units would reveal that g(0, t) does in fact approach a constant value thus indicating a steady displacement. In case of p2 = 0.7 we can make the same conclusion by setting the initial interface to be elongated so that g(0,0) is greater than eventual value to which g(0,t) converges and observe g(0,t) decreases from that value (an idea analogous to that used in lubrication model). Plot of an elongated initial interface and the subsequent track of the interface on the wide side is given in Fig. 5.8. However, for say p2 = 0.8 guessing this eventual value to which g(0,t) converges is not straightforward as well as the dimension of computational domain in the £ direction re-quired would be considerably larger. Evolving interface from an initial horizontal configuration as we have done for smaller values of p2 is also not feasible in this case as the "development time" would be too long. On the other hand, it is clear that for large values of p2 - (1.1 or 1.0) displacement is go-ing to be unsteady with g(0, t) steadily increasing with a large gradient. Thus we can conclude that transition between steady and unsteady displacement according to the 2-D flow simulations occurs somewhere between p = 0.7 and p2 = 1.0 in this case, however acquisition of a more precise bound appears to be problematic. At the very least we can conclude that prediction for 67 Chapter 5. Lubrication model ) 1 1 1 1 1 1 -I 1 . 1 , 1 , H 1 1 1 1 1 1 2 3 4 5 6 1 2 3 4 5 6 2 4 6 8 10 t t t Figure 5.9: Top: plot of the differential velocity WitW — Wi,h calculated using the lubrication model versus the consistency of the displaced fluid k2 with other parameters fixed (e = 0.3; k\ = 1.0; mi = rn.2 = 1-0; Ty\ = 1.0; r y 2 = 0.8; p\ = 1; p2 = 0.95; (3 = 0). On the plot, solid circles indicate the values of consistency for which 2-D flow simulation indicated an unsteady displacement and empty circles where a steady displacement was achieved. Below are plots of the interface position at the wide side of the annulus (4> = 0) from the simulated flow in frame of reference moving with the average speed flow for (left to right): a) k2 = 0.4; b) k2 = 0.6; c) k2 = 1.2. steady displacement based on lubrication model appears to be conservative compared to flow simulation results, i.e. flow simulation predicts a steady displacement for a larger range of p2. Second example illustrating the same kind of comparison is given in Fig. 5.9. Here we fix all parameters at values specified except for the consistency of the displaced fluid k2. Again, similar sort of problems as in the previous example arise with a more exact determination of the transition from steady to unsteady displacement using the flow simulation. Nevertheless we note again that the transition occurs later than predicted by the lubrication model. 68 Chapter 5. Lubrication model 1.0 1.0 0.5" 0.5' -1 , . 1 1 1 1 J 1 1 1 1 1 0 .9 1 .0 1.1 1.2 1.3 1.4 1.0 1.2 1.4 1.6 1.8 p2 k2 Figure 5.10: Plots of the velocity at the wide side of the annulus ((f) = 0) as given by the lubrication model (solid line) and from 2-D flow simulations (crosses). Left to right: a) density of displaced fluid p2 between 0.9 and 1.3 and other parameters fixed at e = 0.4; fci = 1.0; k2 = 0.8; m\ = m2 = 1.0; T y i = 1-0; ry2 = 0.8; p\ = 1; 8 = 0 b) consistency of the displaced fluid k2 between 1.0 and 1.8 with other parameters fixed at e = 0.3; k\ = 1.0; m\ = m2 = 1.0; T y i = 1.0; TV2 = 0.8; p i = 1; p2 = 0.95; 3 = 0. Second characteristic on which we compare these two models is the prediction of the veloc-ity on the wide side of the annulus. It is sensible to make this comparison when we expect an unsteady displacement - as otherwise, in a case of steady displacement, during a simulation an initially elongated interface would level off with velocity on the wide side eventually diminishing to zero, making its estimation ambiguous. We adopt the same strategy as above: for the first example we fix all parameters except for density of the displaced fluid p2. We vary p2 between 0.9 and 1.3 and fix other parameters as before at e .= 0.4; h = 1.0; k2 = 0.8; mx = m2 = 1.0; ry\ = 1.0; T y 2 = 0.8; pi = 1; B = 0. Results obtained thus from the flow simulation together with ones predicted by the lubrication model are given in Fig. 5.10 a). We repeat the procedure by now varying the consistency of the displaced fluid k2 between 1.0 and 1.8 with other parameters fixed at e = 0.3; k\ = 1.0; mi = m2 = 1.0; TYI = 1.0; ry2 = 0.8; p\ = 1; p2 = 0.95; 8 = 0. Results are plotted in Fig. 5.10 b). In both cases we note that the estimation of speed on wide side by the lubrication model is slightly higher than that obtained from the flow simulations. 69 Chapter 6 C o n c l u s i o n a n d D i s c u s s i o n In this thesis we have described and analyzed in detail a model of non-Newtonian displacement flow in annular geometry. Although building extensively on previously done work in this area we have achieved the following main new results of importance: 1. We derived an immiscible fluids model, with appropriate jump conditions. This model is easier to analyze than the concentration equations formulation. The two models lead to an equivalent weak formulation. 2. We obtained an exact solution for the form of the interface and the stream function in the case of a concentric annulus from the immiscible fluids formulation. This provided a useful validation tool for later results as well as serving as the basis for the perturbation solution. 3. We derived a perturbation solution for the form of the interface and the stream function in the case of small eccentricity e << 1 from the immiscible fluids formulation. This allowed us to easily demonstrate several important traits of the interface shape, depending on variations of parameters such as densities, consistencies and yield stresses of the fluids. 4. We showed existence and uniqueness of solution * for the stream function equations for a given concentration field for the concentration formulation (or equivalently a given interface position for the immiscible fluids formulation). 5. We implemented an iterative algorithm for the solution of the concentration formulation of the stream function equations that solves the "exact" problem without making additional 70 Chapter 6. Conclusion and Discussion simplifying assumptions in the model or regularization of the constitutive equations. We examined convergence of the algorithm and validated it against analytic solutions. 6. We derived a lubrication type model from the immiscible fluids formulation that can reliably predict the flow type - differentiating between stable and unstable displacement and formation of a mud channel. The associated computations take milliseconds of C P U time to complete and compare well with those produced by 2-D flow simulation. There are several comments that should be made regarding the material presented here. As mentioned earlier, we take the flow model summarized in [15] as the staring point for our analysis. As such, limitations of that model all apply here - and one should consult [15] for a more comprehensive discussion of these. The primary problems and limitations outlined are i) assumption of homogeneity of fluids across the gap; ii) assumption of non-turbulent flow in the whole annular domain - in practice some cement jobs are pumped in turbulent regime; iii) ignoring via scaling laws what happens close to the interface and concentrating on the bulk flow instead; iv) assumption of an imposed flow rate throughout the well, which can break down if U-tubing occurs; v) assumption of the casing being stationary - sometimes the casing is slowly rotated during the cementing process which has the effect of shear thinning both fluids. Some of these can be clearly addressed in the future by expanding the model while the importance of others can be examined experimentally. In particular examination of assumption (i) can be carried out by extending the current research on displacements in a long axial section of the annulus, considered in [1, 16]. For the augmented Lagrangian method implemented here for the solution of the stream function, the one major limitation appears to be the computational speed. For a 20 x 100 computational domain considered here a typical flow simulation involving about 3000 time steps (equivalent to 30 time units, with time step ~0.01 obtained from C F L condition) can take several hours on current mid-range PCs depending on parameters. This can be remedied by a more extensive exploitation of the over-relaxation parameters pn and r, implementation of more efficient nu-merical methods for solution of the algorithm's equations or using a non-regular mesh to reduce 71 Chapter 6. Conclusion and Discussion the number of control volumes needed. Lastly we note the relation of lubrication model to the 2-D flow simulations. The predic-tions made by the lubrication model are, as we would expect, conservative. This can be seen in Fig. 5.7 and 5.9 where the domain on which the lubrication model predicts a steady flow is in all cases smaller than that on which the 2-D flow simulation predicts a stable displacement. Consistently, as shown Fig.5.10 the velocity of the interface at the wide side of the annulus pre-dicted by the lubrication model is greater than that computed from the 2-D flow simulations. Determination of just how conservative is the lubrication model is complicated by the difficul-ties in determining the flow type from the 2-D simulations described in 5.2.3 using a particular example. Some indication of the order of magnitude of the discrepancy however can be perhaps deduced from comparing the predicted speeds at the wide/narrow sides of the annulus. 72 B i b l i o g r a p h y M. Allouche, LA. Frigaard, and G. Sona. Static wall layers in the displacement of two visco-plastic fluids in a plane channel. J. Fluid Mech., 4%4, pages 243-277, 2000. N.A. Barton, G.L. Archer, and G.L. Seymour. Computational fluid dynamics improves liner cementing operation. Oil and Gas Journal, pages 92-98, September 1994. R.M. Beirute, F.L. Sabins, and K.V. Ravi. Large-scale experiments show proper hole con-ditioning: A critical requirement for successful cementing operations. Society of Petrolium Engineers, 1991. Paper number SPE 22774. S.H. Bittleston. Laminar displacement in inclined eccentric annulus. Confidential Schlum-berger report SCR/SR/1990/003/FLM/C, (1990). S.H. Bittleston. Laminar displacement in vertical eccentric annulus. Confidential Schlum-berger report SCR/SR/1990/032/FLM/C, (1990). S. Brady, P.P. Drecq, K.C. Baker, and D.J . Guillot. Recent technological advances help solve cement placement problems in gulf of mexico. Society of Petrolium Engineers, 1992. Paper number IADC SPE 23927. P. Coussot. Saffman-taylor instablity in yield-stress fluids. J. of Fluid Mech., 1999. M. Couturier, D. Guillot, H . Hendriks, and F. Callet. Design rules and associated spacer properties for optimal mud removal in eccentric annuli. Society of Petrolium Engineers, 1990. Paper number SPE 21594. R.J. Crook, S.R. Keller, and M.A. Wilson. Solutions to problems associated with deviated-wellbore cementing. Society of Petrolium Engineers, 1985. Paper number SPE 14198. J. Ferguson. The placement simulator: laminar displacement module. Confidential Schlum-berger report SCR/SR/1990/039/FLM/C, (1990). R.W. Flumerfelt. An analytical study of laminar non-newtonian displacement. Society of Petrolium Engineers, 1973. Paper number SPE 4486. M. Fortin and R. Glowinski. Augmented Lagrangian Methods. North-Holland, 1983. LA. Frigaard. A quick and robust 2-d simulator for laminar fluid-fluid displacements in primary cementing. Confidential Schlumberger report number 098-016, (1998). LA. Frigaard. Wellbore fluid-fluid displacement research and engeneering review, 1997. LA. Frigaard, S.H. Bittleston, and J. Ferguson. Mud removal and cement placement during primary cementing of an oil well. Submitted to J . Engng. Mathematics 2001, to appear. [16] LA. Frigaard, O. Scherzer, and G. Sona. Uniqueness and non-uniqueness in the steady displacement of two visco-plastic fluids. ZAMM, 81(2), pages 99-118, 2001. 73 Bibliography R. Glowinski. Numerical Methods for Nonlinear Variational Problems. Springer-Verlag, 1983. D. Guillot, D. Hendriks, F. Callet, and B. Vidick. Well Cementing, chapter Mud Removal (Chapter 5). Schlumberger Educational Services, Houston, 1990. R.C. Haut and R.C. Crook. Primary cementing: The noncirculatable mud displacement process. Society of Petrolium Engineers, 1979. Paper number SPE 8253. R.C. Haut and R.J . Crook. Primary cementing: Optimizing for maximum displacement. World Oil, pages 105-116, November 1980. G. M . Homsy. Viscous fingering in porous media. Ann. Rev. Fluid Mech., 1987. R.R. Huilgol and M.P. Panizza. On the determination of the plug flow region in bingham fluids through the application of variational inequalities. J. Non-Newtonian Fluid Mech., 58, 1995. A . Jamot. Displacement de la boue par le latier de ciment dans l'espace annulaire tubage-paroi d'un puits. Revue Assoc. Franc. Techn. Petr., 1974. 224, pp. 27-37. S.R. Keller, R .J . Crook, R.C. Haut, and R.C. Kulakofsky. Problems associated with deviated-wellbore cementing. Society of Petrolium Engineers, 1983. Paper number SPE 11979. V . C . Kellessidis, R. Rafferty, A . Merlo, and R. Maglione. Simulator models u-tubing to improve promary cementing. Oil and Gas Journal, March 1994. A . Lindner, P. Coussot, and D. Bonn. Viscous fingering in a yield stress fluid. Physical Review Letters, 85 Vol. 2, 2000. C.F. Lockyear and A.P. Hibbert. Integrated primary cementing study defines key factors for field success. Journal of Petrolium Technology, December 1989. C.F. Lockyear, D.F. Ryan, and D.F. Gunningham. Cement channeling: How to predict and prevent. Society of Petrolium Engineers, 1989. Paper number SPE 19865. M . Martin, M . Latil, and P. Vetter. Mud displacement by slurry during primary cementing jobs - predicting optimum conditions. Society of Petrolium Engineers, 1978. Paper number SPE 7590. R .H . McLean, C W . Manry, and W.W. Whitaker. Displacement mechanics in primary cementing. Society of Petrolium Engineers, 1966. Paper number SPE 1488. M . Mineev-Weinstein. Selection of the saffman-taylor finger width in the absence of surface tension: An exact result. Physical Review Letters, Vol. 80, pages 2113-2116, 1998. H . Pascal. Dynamics of moving interface in porous media for power law fluids with yield stress. Int. J. Engng. Sci., Vol 22 No.5, pages 577-590, 1984. H . Pascal. Rheological behaviour effect of non-newtonian fluids on dynamic of moving interface in porous media. Int. J. Engng. Sci., Vol 22 No.3, 1984. 74 Bibliography [34] H . Pascal. Rheological effects of non-newtonian behaviour of displacing fluids on stability of a moving interface in radial oil displacement mechanism in porous media. Int. J. Engng. Sci., Vol 24 No.9, pages 1465-1476, 1986. [35] H . Pascal. A theoretical analysis on stability of a moving interface in a porous medium for bingham displacing fluids and its application in oil displacement mechanism. The Canadian Journal of Chem. Engng., Vol 64, pages 375-379, 1986. [36] K . V . Ravi, K . V . Beirute, and R.L. Covington. Erodability of partially dehydrated gelled drilling fluid and filter cake. Society of Petrolium Engineers, 1992. Paper number SPE 24571. [37] K . V . Ravi, R . M . Beirute, and R.L . Covington. Improve primary cementing by continious monitoring of circulatable hole. Society of Petrolium Engineers, 1993. Paper number SPE 26574. [38] D.F. Ryan, D.S. Kellingray, and C F . Lockyear. Improved cement placement on north sea wells using a cement placement simulator. Society of Petrolium Engineers, 1992. Paper number SPE 24977. [39] P.G. Saffman. Viscous fingering in hele-shaw cells. J. Fluid Mech., Vol. 173, pages 73-94, 1986. [40] P.G. Saffman and G. Taylor. The penetration of a fluid into a porous medium or hele-shaw cell containing a more viscous fluid. Proc. Roy. Soc. London, Vol. 245, pages 312-329, 1958. [41] C.W. Sauer. Mud displacement during cementing: A state of the art. Journal of Petrolium Technology, pages 1091-1101, September 1987. [42] A . Tehrani, S.H. Bittleston, and P.G.J . Long. Flow instabilities during annular displace-ment of one non-newtonian fluid by another. Experiments in Fluids 14, pages 246-256, 1993. [43] A . Tehrani, J . Ferguson, and S.H. Bittleston. Laminar displacement in annuli: A combined theoretical and experimental study. Society of Petrolium Engineers, 1992. Paper number SPE 24569. [44] E .H . Vefring, K.S. Bjorkevoll, S.A. Hansen, N . Sterri, O. Saevareid, 0 . Aas, and 0 . Merlo. Optimization of displacement efficiency during primary cementing. Society of Petrolium Engineers, 1997. Paper number SPE 39009. 75 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080071/manifest

Comment

Related Items