TRANSIENT EFFECTS IN OILFIELD C E M E N T I N G FLOWS by MIGUEL ANGEL MOYERS GONZALEZ B.Sc, Instituto Tecnologico Autonomo de Mexico, 2000 M.Sc:, The University of British Columbia, 2002 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in THE FACULTY OF GRADUATE STUDIES (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA January 2006 © Miguel Angel Moyers Gonzalez, 2006 Abstract In this thesis we study instabilities and mud channel removal during the primary cementing of an oil well. This process involves displacement of a sequence of non-Newtonian fluids along a narrow eccentric annulus. Previously, this has been modelled as a pseudo-steady process using a Hele-Shaw approximation. In such models, the stream function is governed by a steady elliptic problem and the fluids, (modelled via a concentration equation), simply advect along the annulus. It has been shown that for certain rheological and physical parameters, the displacement front will advances much faster on the wide side of the annulus than on the narrow side. In extreme cases the displacement front does not advance at all on the narrow side of the annulus and a static mud channel results as thefingeradvances up the wide side. In this thesis we consider whether the interface of a progressively advancing finger will remain stable to small perturbations. There is in fact experimental evidence that interfacial instabilities can occur in this situation. We find that the interface is in fact stable whenever there is a static mud channel on the narrow side of the annulus. Consequently, we also investigate how a mud channel might be removed by pulsation of the flow rate. Study of these two phenomena cannot be undertaken with the pseudo-steady framework. Therefore, we extend this model to flows that are fully transient. The transient model consists of a nonlinear evolution equation for the stream function. In chapter 3 we show that this transient model is in fact well-posed. In chapter 4 we study stability of multi-layer parallel flows, i.e. long fingers. If both fluids are yielded at the interface, instabilities may arise for different combinations of the 14 dimensionless parameters. These instabilities are related to a jump in tangential velocity at the interface and do not appear to have been identified before. In chapter 5 we investigate the case where a static mud channel develops on the narrow side of the annulus. Our stability theory predicts only linear stability. We therefore study the effects of a finite pulsation of the flow rates via numerical simulation. It seems that if we perturb the flow from the beginning of the displacement, the transient model fully captures the effects of the perturbation and the width of the mud channel is reduced. The pseudo-steady velocity model does not report any significant changes with respect to the results using a constant flow rate. If however, pulsation is applied after the mud channel has already formed, the removal of the mud channel will be unsuccessful. ii Table of Contents Abstract ii Table of Contents iii List of Tables vi List of Figures vii Acknowledgements xii xiii Chapter 1. Introduction 1 1.1 Primary cementing process description 1.2 Non-Newtonian Fluids 1.2.1 Time independent non-Newtonian 1.3 Literature Review 1.3.1 Industrial background 1.3.2 Displacement flow in different geometries 1.3.3 Displacements in a Hele-Shaw cell 1.3.4 Interfacial instabilities 1.4 Overview Chapter 2. 1 3 \ 3 5 5 7 9 13 14 fluids Model derivation 16 2.1 The pseudo-steady model and some open questions 2.1.1 Open problems and limitations 2.2 Modelling transient bulk flow cementing displacements 2.2.1 Non-Dimensionalisation 2.2.2 Reduced shear flow equations 2.2.3 2D averaged velocity field 2.2.4 Fluid concentration field 2.3 Reduced shearflowmodels 2.3.1 Stability of u and u , at (</>,£): 2.3.2 Stability of u and u , at (</>, £): 2.4 Eccentric annular Hele-Shaw displacement transient model 2.4.1 Rheological assumptions 2.4.2 Hele-Shaw model derivation 2.4.3 Boundary conditions 2.5 Interface tracking model derivation 17 19 20 22 24 25 26 29 31 33 35 35 37 40 42 ; s a Chapter 3. s Qualitative results . 3.1 Existence & uniqueness of ^ 3.1.1 Preliminary results 47 48 49 iii Table of Contents 3.1.2 3.2 3.3 Definition of dI 2[u] hL 3.1.3 Steady state problem 3.1.4 Transient problem 3.1.5 Existence of a solution and dI ^M 56 • • • • 57 59 Continuity 60 3.2.1 Preliminary results 60 69 3.2.2 Continuity with respect to rheological and physical properties 3.2.3 Steady State 70 3.2.4 Continuity for Evolution equation 76 Qualitative behaviour of solutions 79 3.3.1 E n d conditions and test concentrations 80 3.3.2 Steady state problem 81 3.3.3 Transient problem 84 Chapter 4. Stability of multi layer 4.1 4.2 55 2 flows 90 Parallel multi-layer flows 92 4.1.1 Static mud channels 97 4.1.2 Summary of results 101 Stability of parallel multi-layer flows 101 4.3 Single fluid problem 103 4.4 Two-fluid eigenvalue problems 106 4.5 Analytical simplifications 110 4.5.1 Newtonian fluids in a concentric annulus 110 4.5.2 Power-Law fluids in a concentric annulus 4.6 119 4.6.1 119 . 4.6.2 4.7 4.8 Numerical Method Yield stress fluids 123 Long wavelength asymptotics 127 4.7.1 Newtonian fluids 127 4.7.2 Yield stress fluids 131 4.7.3 Industrial application 138 Spatial stability problem 141 4.8.1 Single fluid spatial instability 142 4.8.2 Two fluids spatia! instability 142 Chapter 5. Static mud channels 5.1 5.2 115 Numerical Results 145 Augmented Lagrangian Method 146 5.1.1 Principles of the method 146 5.1.2 ALG2 148 5.1.3 Convergence of A L G 2 149 Application of A L G 2 to the displacement problem 151 5.2.1 Steady model 151 5.2.2 Transient model 154 Finite Element discretization 155 5.4 Flux Corrected Transport Algorithm 156 5.5 Numerical simulations for the annular displacement with pulsating flow rate 160 5.3 iv Table of Contents 5.5.1 Validation 5.5.2 Pulsation effects Chapter 6. ...160 161 Summary and conclusions 183 6.1 Main contributions of the thesis 6.2 Interfacial instabilities 6.3 Mud channel removal 183 184 185 Bibliography 187 v List of Tables 4.1 Code validation List of Figures 1.1 1.2 Schematic of the primary cementing process: a) Drilling to desired depth, b) Placement of the steel casing, well bore still full of drilling mud. c) Cement is pumped down inside of the steel casing, d) Cement displaces drilling mud through the annular region, e) Cement hardens into hydraulic seal for zone isolation 2 Schematic behaviour of shear stress vs. shear rate for time independent non-Newtonian fluids 4 2.1 2.2 Geometry of the well Example of concentration field at fixed £ and t 3.1 The function x(|V*|,f/',7v,K;,// ,n) 51 a) 7 > T T / 2 , b) (3< T T / 4 , 7 < T T / 2 , C ) R > T T / 4 , 7 < T T / 2 .0 < T T / 4 , 7 < T T / 2 65 a) Decay of ||* — ^s\\L (n)- b) Fixed concentration field 88 Transient decay to steady state solution, contour plots of \I> at time t = 1, 5 , 10, and contour plot of the steady state solution \I/ 89 3.2 3.3 3.4 3.5 18 45 < 00 6 4 2 S 4.1 4.2 4.3 4.4 4.5 4.6 Multilayer baseflow,both fluids yielded at the interface. Rhelogical and physical parameters are: 7V, 1 = 2, K\ = 2, p\ = 2, mi = 2, 7V,2 = 1, «2 = 1, P2 = 1, mi = 1, e = 0.4, fa = 0 . 5 and P = n/W '. 95 Multilayer base flow, bothfluidsyielded at the interface, mud channel in the narrow side of the annulus. Rhelogical and physical parameters are: 7v,i = 0.5, K\ = 1, p\ = 2 , mi = 2 , 7v,2 = 0 . 6 , «2 = 1, pi = 1, m2 = 1, e = 0 . 5 , fa = 0 . 5 and P = TT/10 95 Multilayer baseflow,fluid2 unyielded at the interface. Rhelogical and physical parameters are: 7v,i = 1, KI = 1, pi = 4, mi = 2 , TV,2 = 10, K2 = 1, P2 = 1, rn,2 = 1, e = 0.8, fa = 0 . 5 and P = T T / 1 0 96 Multilayer baseflow,fluidl unyielded at the interface. Rhelogical and physical parameters are: 7v,i = 7, K\ = 7, p\ = 1, mi = 2, 7v,2 = 0.2, «2 = 1, P2 = 1, m2 = l, e — 0.2, fa = 0 . 5 and p = 0 96 Contour plot of fa in, eccentricity vs. yield stress offluidone, buoyancy parameter equal zero. Rhelogical and physical parameters are: K\ — 0.5, pi = 1, mi — 1, 7v,2 = 0.5, K2 = 0.4, P2 = 1, mi = 1 and P = 0. The zero level curve corresponds to the equality of condition ( 4 . 1 4 ) , i.e. fa mm — 1) the rest of the level curves correspond to different values of ^>i,min when condition ( 4 . 1 4 ) is fulfilled 99 Contour plot of fa, ini eccentricity vs. yield stress offluidone, buoyancy parameter negative. Rhelogical and physical parameters are: «i = 0.5, pi = 1.1, mi = 1, 7v, = 0.5, = 0.4, pi — 1, mi = 1 and p-= 0. The zero level curve corresponds to the equality of condition ( 4 . 1 4 ) , i.e. fa i — 1, the rest of the level curves correspond to different values of </>i,min when condition ( 4 . 1 4 ) is fulfilled. 99 >m m 2 K 2 im n vn List of Figures 4.7 Contour plot of <f>i,min, eccentricity vs. yield stress of fluid one, buoyancy parameter positive. Rhelogical and physical parameters are: K\ — 0.5, p\ = 1, mi = 1, 7y,2 = 0.5, « 2 = 0.4, p = 1.1, m = 1 and 0 — 0. The zero level curve corresponds to the equality of condition (4.14), i.e. (p^ m = 1, the rest of the level curves correspond to different values of 0i,min when condition (4.14) is fulfilled 100 Contour plot of 4>i,wini eccentricity vs. yield stress of fluid one, equal consistency. Rhelogical and physical parameters are: K\ = 0.5, pi = 1, mi = 1, 7y — 0.5, K = 0.5, P2 = 1, m = 1 and 0 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. 0i, in = 1, the rest of the level curves correspond to different values of <&i, i when condition (4.14) is fulfilled '. 100 Contour plot of 4>^ im eccentricity vs. yield stress of fluid one, consistency of fluid 1 less than consistency of fluid 2. Rhelogical and physical parameters are: K\ = 0.3, pi = 1, mi = 1, TY,2 = 0.5, K2 = 0.4, p = 1, m = 1 ^ 0 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. <fo, j = 1, the rest of the level curves correspond to different values of 0j i when condition (4.14) is fulfilled 101 Spectrum of the temporal normal mode with lower and upper bounds as in (4.31) & (4.32), figures on the left-hand-side show a fluid with rheological parameters: left-handside, 7y,i = 1, KI = 1, pi = 0.5, mi = 1, e = 0.4 and 0 = 0 on the right-hand-side we only change pi to 1 105 Maximum imaginary part of s vs. a. Rheological and physical parameters are: m = 1, e = 0.0, (pi = 0.5 and 0 = 0. a) pi = 1, p = 0.75. b) pi = 1, p = 1.25 112 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1, = l, = 1 = 0.0, & = 0.5 and 0 = 0 112 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1, e = 0.0, <fc = 0.45 and /? = 0. a) p = 1, p = 0.75. b) p = 1, p = 1.25 113 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI — 1, e = 0.0, & = 0.55 and 0 = 0. &) pi = 1, p = 0.75. b) pi = 1, p = 1-25 113 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1, e = 0.0, (pi = 0.45 and 0 = TT/8. a) pi = 1, p = 0.75. b) pi = 1, p = 1.25 114 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1.5, K = 1, e = 0.0 and /? = 0, mi = 1.5, m = 2. a) = 0.5. b) <pi = 0.6. For b < 0, pi = 1 and p = 0.75,forb = 0, pi = 1 and p = 1 andforb > 0, pi = 1 and p = 1.25. .. 117 Maximum imaginary part of s vs. a. Rheological and physical parameters are: Ki = 1.5, K2 = 1, e = 0.0, 0i = 0.5 and /? = TT/8, mi = 1.5, m = 2. For 6 < 0, pi = 1 and p2 = 0.75, for b = 0, pi = 1 and p = 1 andfor6 > 0, pi = 1 and p = 1.25. . . . . . . . 118 Maximum imaginary part of s vs. a. Inclined annulus. Rheological and physical parameters are: KI = 1.5, K = 1, mi = 1, m = 1, b = —1 and 0 = 0. a) Exact solution, b) Numerical solution 121 Maximum imaginary part of s vs. a. Rheological and physical parameters are: «i = 1.5, « 2 = 1, e = 0.0 d 0 = 0, mi = 1, m = 1, ry^ = 0.9 and 7y, = 0.7. a) (pi = 0.5. b) (pi = 0.6. For 6 < 0, pi = 1 and p = 0.75, for b = 0, pi = 1 and p = 1 and for b > 0, pi = 1 and p = 1.25 123 2 2 m 4.8 ]2 2 2 m 4.9 m m a n 2 2 mn im n 4.10 4.11 2 4.12 P l 4.13 P 2 e 2 x 4.14 2 x 2 4.15 2 2 2 4.16 2 2 2 4.17 2 2 2 2 2 4.18 2 2 4.19 2 a n 2 2 2 2 2 viii n List of Figures 4.20 Maximum imaginary part of s vs. a. Rheological and physical parameters are: a) KI = 1.5, K2 = 1, e = 0.0, fa = 0.5 and (3 = 0, mi = 1, m = 1, TY,I = 0.7 and 7V,2 = 0.9. b) K\ — 1.5, K = 1, = 0.0, fa = 0.5, 3 = IT/8, mi = 1, m = 1, 7v,i = 0.7 and 7 v , 2 = 0.9., mi = 2, m,2 = 1-5. For b < 0, pi = 1 and p = 0.75, for 6 = 0, p\ = 1 and P2 = 1 and for 6 > 0, pi = 1 and p = 1.25 124 4.21 Maximum imaginary part of s vs. a. Different positions of the interface between the fluids. Rheological and physical parameters are: K\ — 1.5, K = 1> e = 0.0 and /3 = 0, mi = 1, m = 1, TV,I = 0.9 and TV,2 = 0.7. a) e = 0.2. b) e = 0.4. c) e = 0.6. d) e = 0.8. For b < 0, pi = 1 and p = 0.75, for b = 0, pi = 1 and p = 1 and for b > 0, pi = 1 and p = 1.25 125 4.22 Spectrum and eigenfunctions for system (4.39)-(4.40), a = 0.5. Rheological and physical parameters are: K\ = 1.5, K2 = 1, mi = 1, m = 1, TV,I = 0.9, 7v,2 = 0.7, negative buoyancy, e = 0.4 and j3 — 0. a) Spectrum of the eigenvalue problem, b) Eigenfunctions of the interfacial mode, solid line shows the real part of fe and / , dashed line shows the imaginary part of /1 and fe and dashed thick line shows the magnitude of /1 and fe. c) Eigenfunctions of an interior mode, solid line shows the real part of /1 and / , dashed line shows the imaginary part offeandfeand dashed thick line shows the magnitude of fi and fe 126 4.23 Contour plot of s . Rheological and physical parameters are: K\ = 1, pi = 1, e = 0.0, fa = 0.5 and P = 0 131 4.24 Contour plot of s , «i vs. K . Rheological and physical parameters are: pi = 1, mi = 1, m = 1, 7v,i — 1, TV,2 = 1» P2 = 1, e = 0.2,fa= 0.5 and P = 0 136 4.25 Contour plot of s , K\ VS. K - Rheological and physical parameters are: mi = 1, m = 1> 7v,i = 1, 7 v , 2 = 1, e = 0.2,fa= 0.5 and P = 0. a)pi = 1.1 and p = 1. b)pi = 1 and p '= 1.1 137 4.26 Contour plot of s , n\ VS. K . Rheological and physical parameters are: pi = 1, TV,I = 1, 7 v , 2 = 1, P2 = 1-1, e = 0.2, (pi = 0.5 and P = 0. a) mi = 1.2 and m = 1. b) mi = 1 and m = 1.2, 137 4.27 Contour plot of s , K,\ VS. K , buoyancy parameter equal zero. Rheological and physical parameters are: KI = 1, pi = 1, /t = 1, p = 1, e = 0.2,fa= 0.5 and P = 0. a) TV,I = 1.5 and TV,2 = 1- b) TV,I = 1 and 7v,2 = 1-5 138 4.28 Example unsteady displacement and local instability/viscous fingering regimes: effects of changing TV,2 and e for different K : (a) « — 0.5; (b) re = 1; (c) K = 1.5. Fixed parameters are: TV,I = 1-0, K.\ — 1.0, mi = 1, m = 2, b = —0.5, P = 0. F, local fingering instability ; U, unsteady, but no local instability; S, neither unsteady. From Pelipenko & Frigaard (2004) 140 4.29 Maximum imaginary part of s vs. a. Rheological and physical parameters are: n\ = 1, 7v,i = 1, mi = 1, m = 2, b = -0.5, fa = 0.5 and P = 0. a) K = 0.5. b) K — 1- c) K = 1-5 141 2 e 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 5.1 Convergence of computational solution to stable steady state displacement: a) Successive contours of c(faz,t) = 0.5 at times t = 0, 0.2, 0.4, ... 6. b) Analytical solution for the steady displacement, c) Convergence of the wide and narrow side interface positions to the steady state. Rheological and physical parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, pi = 1, p = 0.9, TV,I = 1, TV.i = 0.9, TV,I = 0.7 e = 0.0 and P = 0 161 Initial condition for the displacements 162 2 2 5.2 2 ix List of Figures 5.3 Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2-K/UI, U = 10 and S = 0. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, 7V,i = 1, T V , I = 0.9, 7y,i = 0.7 e = 0.3, 0 = 0 165 5.4 Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2ix/u), u> = 10 and S = 0.1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: n\ = 0.5, K = 0.4, mi = 1, m = 1.2, pi = 1, p = 0.9, TV,I = 1, T V , I = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 166 5.5 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2-ir/wt, ui = 10, 5 = 0.1 and e = 0.1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, pi = 1, p = 0.9, ry,i = 1, T V , I = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 167 5.6 Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2-K/LOC, to = 10, S = 0.1 and e = 1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, ry,i = 1, T V , I = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 168 5.7 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2-K/IJOC, W = 10, 5 = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: KI = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, P2 = 0.9, T V , I = 1, iY,i = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 169 5.8 Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2w/uje, w = 1, S = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: n\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, ry,i = 1, T V , I = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 170 5.9 Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 27r/we, UJ = 0.1, 5 = 0.1 and e = 10. Interface position and velocity fieldfortimes: a)-b) T/4. Physical and rheological parameters: KI = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, 7v,i = 1, T V , I = 0.9, 7v,i = 0.7 e = 0.3, 0 = 0 171 5.10 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2ir/u>, u = 10 and 5 = 0. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: «i = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, 7v,i = 1, T V , I = 0.9, 7v,i = 0.7 e = 0.5, 0 = 0 172 5.11 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2TT/LJ, UJ = 10 and 5 = 0.1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: n\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, 7v,i = 1, T V , I = 0.9, 7Y,i = 0.7 e = 0.5, 0 = 0 173 p 2 2 2 p 2 2 2 P 2 2 2 2 2 p 2 P 2 2 p 2 2 2 p 2 2 2 P 2 2 2 P 2 2 x 2 List of Figures 5.12 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2rr/W, w = 10, 5 = 0.1 and e = 0.6. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: « i — 0.5, K = 0.4, mi = 1, m — 1.2, p i = 1, p = 0.9, 7 V , i = 1, T V , ! = 0.9, T V , I = 0.7 e = 0.5, 8 = 0 ' 174 , 5.13 Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Steady state velocity model. Period T = 2TT/LJ, W = 10 and 8 = 0. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) -3T/4. g)-h) T. Physical and rheological parameters: KI = 0.5, K = 0.4, mi = 1, m = 1.2, p i = 1, p = 0.9, 7v,i = 1,' T V , ! = 0.9, T V , I = 0.7 e = 0.8, 8 = 0 175 5.14 Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Steady state velocity model. Period T = 2TT/W, LJ = 10 and S — 0.1. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p i = 1, p = 0.9, 7v,i = 1, T V , I = 0.9, 7 v , i = 0.7 e = 0.8, 8 = 0. 176 5.15 Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Steady state velocity model. Period T = 2-K/LJ, LJ = 10 and 5 = 0.2. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: « i = 0.5, K = 0.4, mi = 1, m = 1.2, p i = 1, P2 = 0.9, TV,I = 1, T V , I = 0.9, T V , I = 0.7 e = 0.8, 8 = 0 177 5.16 Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Transient velocity model. Period T = 2ir/ue, u = 10, 6 = 0.1 and e = 0.6. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ — 0.5, K2 = 0.4, mi = 1, m = 1.2, p i = 1, p = 0.9, T V , ! = 1, 7v,i = 0.9, T V , I = 0.7 e = 0.8, 8 = 0 178 5.17 Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Transient velocity model. Period T = 2ir/ue, UJ = 10, 5 = 0.2 and e = 0.6. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, « — 0.4, mi = 1, m = 1.2, p i = 1, p = 0.9, T V , I = 1, T V , I = 0.9, T V , I = 0.7 e = 0.8, B = 0 179 5.18 Displacement flow in eccentric annulus, interface propagation. Mud channel. Steady state velocity model. Period T = 2W/LJ, LJ = 10 and 5 = 0.2. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, /t = 0.4, mi = 1, m = 1.2, p i = 1, p = 0.9, TV,I = 1, TV,I = 0.9, T V , I = 0.7 e = 0.8, (3 = 0 181 5.19 Displacement flow in eccentric annulus, interface propagation. Mud channel. Transient velocity model. Period T = 2n/u>e, LJ = 10, 5 = 0.2 and e = 0.6. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p i = 1, p = 0.9, TV,I = 1, T v i = 0.9, T V I = 0.7 e = 0.8, 8 = 0 ' 182 P 2 2 2 P 2 2 2 p 2 2 2 p 2 2 ; P 2 2 2 2 p 2 p 2 2 2 p 2 2 xi 2 Acknowledgements I would like to thank Dr. Ian Frigaard, my supervisor, for all the support, guidance and enormous patience he has showed during my stay at UBC. But most of all, for his friendship and trust, to have those I am honored. Certainly, without him this thesis would not be completed. I wish to thank as well, the members of my research committee, Dr. Neil Balmforth, Dr. James Feng and Dr. Tai-Peng Tsai for their useful advice and comments. Dr. Otmar Scherzer deserves my full gratitude for a great and fruitful stay at University of Innsbruck. Also I would like to thank my family and friends back home, although we are many miles away their support and good wishes were with me every day. My family here in Vancouver is one of the main reasons why I am able to finish this period of my life. Rodrigo, Greta, Rafa and Mariana, without you I would not be here now, thank you so much. Thanks to all my friends and colleagues at UBC and Vancouver, my stay here has been as pleasant as it could be in good part because of you. To La Morena, for being with me every single moment of my life. This thesis has been funded by NSERC and by Schulumberger Oilfield Services via the collaborative research grant number 245434. I am grateful for this financial support. During the first year of the PhD thesis, I was recipient of a scholarship from CONACYT. I am grateful for the financial support. Lastly, but certainly not the least, I wish to thank my wife, Dani. Without all her support, care and love I would not have been able to go through the last months of writing and tidying up the thesis. I do not have words to thank her enough. xii To my parents and my wife xiii Chapter 1 Introduction The aim of this thesis it to study instabilities and mud channel removal during the primary cementing of an oil well. This process involves displacement of a sequence of non-Newtonian fluids along a narrow eccentric annulus. It has been shown that for certain rheological and physical parameters of the fluids, the displacement front will advance much faster on the wide side of the annulus than on the narrow side. In extreme cases the displacement front does not advance at all on the narrow side of the annulus and a static mud channel results as the finger advances up the wide side. We consider whether the interface of a progressively advancing finger will remain stable to small perturbations and also investigate how a mud channel might be removed by pulsation of the flow rate, which is sometimes used in the cementing industry. 1.1 Primary cementing process description The primary cementing process proceeds as follows: After drilling the well to the desired depth, Figure 1.1a, the drillpipe is removed and a large steel casing is placed down to the bottom of the well, leaving a gap of « 2cm between the outside of the tube and the inside of the wellbore, i.e. an annulus. At this time, the drilling mud is still in the wellbore, Figure 1.1b. To remove the drilling mud, the cement is first pumped down the inside of the steel casing, Figure 1.1c, and emerges at the end where it then flows into the narrow annular gap between the casing and the well bore, displacing the drilling mud, Figure l.ld. The cement later hardens into an hydraulic seal, Figure l.le. Without the hydraulic seal, well productivity can be reduced significantly and fluids can also leak up to the surface, where they may become an environmental or safety hazard. The fluids used in cementing, (drilling muds, spacerfluidsand cement slurries), are 1 Chapter 1. Introduction Figure 1.1: Schematic of the primary cementing process: a) Drilling to desired depth, b) Placement of the steel casing, well bore still full of drilling mud. c) Cement is pumped down inside of the steel casing, d) Cement displaces drilling mud through the annular region, e) Cement hardens into hydraulic seal for zone isolation . generally characterised as viscoplastic shear-thinning fluids. In practice a sequence of fluids is pumped along an annulus (spacer fluids and cement slurries), each fluid displacing the one in front, and the annular geometry changes slowly in the axial direction. However, the fluid volumes pumped are relatively large and the essential dynamics of the displacement may be studied by considering what happens between any 2 fluids on an annular section of constant geometry. Note that Figure 1.1 presents an ideal case. Usually the angle of inclination of the well is different from zero. Centralizers are attached to the steel casing to prevent the casing from slumping to the lower side of the annulus. Even though they prevent the slumping, there are not always enough centralizers fitted and/or the stiffness of the centralizers is not designed well. Thus, the annulus nearly always is eccentric. Another problem is that a full removal of the drilling mud is not always achieved. A mud channel is frequently left along the narrow part of the annulus. This dehydrates during setting of the cement and the dried mud forms a porous conduit along the length of the well. Therefore, complete zonal isolation is not achieved. 2 Chapter 1. Introduction 1.2 Non-Newtonian Fluids As stated in the previous section, the fluids used in cementing are generally characterized as visco-plastic shear thinning fluids. In this section we give a short summary of the constitutive equations which characterize these fluids. For a more complete review see [7, 15, 37]. Let us assume shear flow of a liquid between two parallel plates which its velocity vector can be expressed as u = (7?/, 0,0) where 7 is the shear rate. In the case of a Newtonian fluid the stress distribution for shear flow is given in the following form: °~yx=' i n cr =a xz yz = 0 a x x -o y y = 0 o- y y -o z z =Q. (1.1) A fluid is called Newtonian if satisfies (1.1) and if rj does not vary with shear rate and is constant with respect to the time of shearing and the stress in the liquid falls to zero after shearing is stopped. A n y liquid that does not satisfies any of these requirements is called non-Newtonian. We will concentrate in fluids for which the difference in the normal stresses is zero, i.e. inelastic fluids, and for which the shear viscosity depends on the shear rate but not explicitly on time. 1.2.1 Time independent non-Newtonian fluids We can subdivide this type of fluids into 3 groups, • Shear thinning (or pseudo-plastic) fluids, for which the viscosity decreases with increasing shear rate. • Shear thickening (or dilatant) fluids, for which the viscosity increases with increasing shear rate. • Yield stress (or visco-plastic) fluids, for which a finite shear stress is required to initiate the flow. 3 Chapter 1. Introduction These fluids are illustrated in Figure 1.2, and we now give some examples of the common fluid models. Shear Rate Figure 1.2: Schematic behaviour of shear stress vs. shear rate for time independent nonNewtonian fluids. For many fluids, as the shear rate approaches zero or infinity the value of 77 approaches constant values, say 770 and T ^ , i.e. these are 2 Newtonian limits. Two models that contain this behaviour, are: V-Voo Vo -Vco (1 + [Ki) ) m (1.2) which is called the Cross model, and V-Voo 1 r?o -Voo (1 + (tf ) ) 2 m/2 (1.3) 7 which is called the Carreau model. Here K and m are empirically fitted constant rheological parameters. If we assume 77 - C 770 and 77 > 7700, at leading order and by a simple redefinition of parameters, the Cross model reduces to 77 = Kj',71-1 4 (1.4) Chapter 1. Introduction which is the well known power law model. If instead we assume that T] -C then again, by a redefinition of parameters, the Cross model reduces to r, = Voo + K i ^ 1 (1.5) which is called Sisko model. If n = 0 in Sisko's model and we redefine the parameters we obtain the Bingham model, (1.6) ri = p + j where p is called the plastic viscosity and ly is the yield stress. The Bingham model is the simplest model that describes visco-plastic behaviour. Other widely used visco-plastic models are: VT = /7y"+ s/K^ v (1.7) which is called the Casson model, and is given in terms of the shear stress, and finally r, = KT~ l + ^ 7 (1-8) which is called the Herschel-Bulkley model. Note that if n = 1 and Ty = 0 we recover the Newtonian model, if n = 1 and iy / Owe recover the Bingham model, and if n < 1 and ry = 0 we recover the power law model. The Herschel-Bulkley model will be the one used for the rest of the thesis. 1.3 1.3.1 Literature Review Industrial background There is an extensive technical literature on primary cementing, see [21] for a review, as this is a large industrial business. Only a part of this literature will be reviewed in this section. Early development of cementing techniques dates back to 1911, we refer to the reader to [65], which gives a good review of the techniques used up to the 1980's. Next we will give a small overview of these techniques. In the 1940's a good amount of attention was given to the cementing process. Investigations on the effectiveness of water pre-flushes ahead of cement, pipe movement, well centered casing, 5 Chapter 1. Introduction mud thinning prior to cementing and turbulent flows were carried out. In the 1960's mud displacements studies with cement slurries in turbulent flow proved to be very effective, this gave to rise to the development of the "scavenger slurry" which is a low density, highly dispersed cement slurry which is pumped before the tail slurry, i.e. the cement that is expected to provide casing support and a hydraulic seal. In the 1970's researchers started looking at how to get a good cementing job in deviated holes. Studies on the effects of centralizers and pipe movement were carried out. In the 1980's attention was given to the design of a mud-cement compatible spacer systems that were capable of transitioning to turbulence at reasonable pump rates. Since the late 60's analytical and experimental studies have been carried out in order to understand better when a good displacement takes place. For example, McLean et.al., [52], considered an analytical model which was later tested against experimental data. They considered the division of the annulus into sectors each one treated as an equivalent sector of a concentric annulus and assumed that the pressure gradient was identical in each sector, (effectively a crude early finite element method!). From this approach they found a critical value for the yield stress of the displacing fluid, in this case cement, for which the pressure gradient in each sector was large enough to initiate flow in all the annular region. This prediction seems to over estimate the actual experimental value of the yield stress of the cement. They claimed that this is due to the fact of neglecting the effect of interfacial shear stresses at the cement mud interface. Later, Clark & Carter, [13] studied the effect of interfacial stresses. They found that these effects can be very important in initiating mud flow and claimed that adverse buoyancy effects could be compensated by displacing the cement at higher rates. Another example of an early analytical model can be found in [9], where the authors approximated the displacement zone as a narrow gap annulus between concentric cylinders, or equivalently the region between two parallel plates. Their model is 1-D pseudo-steady, i.e. they calculate the 1-D steady state solution for the velocity at fixed time. Through integration of the velocity along the interface and determining its position, they are able tofindthe volume of fluid displaced as a function of time. On the other hand, Courtier et.al., [14] actually considered the effect of eccentricity using the basic slot model, which consists in considering independent 6 Chapter 1. Introduction rectangular slots of varying heights at different azimuthal positions. All the analytical models were tested against available experimental data. For the experimental results we refer to the reader to [11, 50, 69, 70, 93]. There have also been a large number of field-based case studies, which we do not review here. For a good displacement in the annular region to happen, almost all laboratory and field studies agree that: • The displacing fluid should be heavier than the displaced one for near vertical annulus. • The use of centralizers to reduce the eccentricity of the annulus because one of the main causes of mud channeling is the casing eccentricity. • Fluids with higher viscosity and yield stress are better able to displace fluids with lower viscosity and yield stress. Most industrial literature agrees that if turbulence is achieved in the annulus it's the best tool to remove the mud from the well. Unfortunately, this is not always attainable in the field. For example the flow rates required may exceed capacity of the pumps at the rig site. 1.3.2 Displacement flow in different geometries We present a brief review of some of the work done in simulating, either analytically, numerically or in a lab setting, the displacement of one miscible fluid by another in different geometries. In 1993, Tehrani et.al., [78], presented experimental results on the displacement of one fluid by another in an eccentric annulus, using both Newtonian and non-Newtonian fluids. Their study shows that the interface between the the two fluids with similar rheologies but different densities is stable as long as the denser fluid lies below the interface. If the interface becomes close to vertical, a bad displacement and a long finger develops. In this case, small perturbations in the flow field may trigger gravity driven instabilities. Azimuthal instabilities also can occur. The intensity of these instabilities is related to density difference and flow rate. These interfacial instabilities are one of the main subjects of the thesis and will be investigated in chapter 4. These results lead to the work reported by Szabo & Hassager in [73]. The fluids they consider 7 Chapter 1. Introduction are Newtonian, their simulations support earlier experimental findings. They solved the full 3-D problem using a mixed Lagrangian-Eulerian method. They solve the equations of motion plus an equation for the moving grid to track the interface between the fluids. Some experimental work concerning the downward or upward vertical displacements of two miscible Newtonian fluids in a vertical pipe can be found in [5, 67]. Scoffoni et.al., [67] considered only the downward displacement of a heavier and more viscous fluid by a lighter and less viscous fluid. A long finger develops along the length of the pipe and at high rates the finger destabilizes showing either an axisymmetric pattern or an asymetric corkscrew mode. Baasubramaniam et.al., [5] considered the opposite case, i.e. the displacement of lighter and less viscous fluids by a heavier and more viscous fluids in downward or upward displacements. For the case of downward displacements they also observed interfacial instabilities. They found that a long finger develops along the length of the duct. It has been shown, [40] that the interface between the fluids is destabilized and thus fingering occurs when a more viscous fluid is displaced by a less viscous one. This is due to the advection of the initial flat interface by the nonuniform velocity profile within the pipe. For downward displacements they also observed interfacial instabilities. For the upward displacements the interface was stable although a "spike" formed at the tip of the finger. Eccentric Hele-Shaw model for annular displacements Bittleston et.al., [10], developed an eccentric Hele-Shaw model for the displacement of one fluid by another in an annular geometry, an oil well. Thefluidsconsidered satisfy the HerschelBulkley model. Because the annulus geometry is long and thin they derive a 2-D model of the bulkfluidmotions by averaging across the gap and scaling the Navier-Stokes equations with the annular gap width, keeping only leading order terms. The model consists of an elliptic partial differential equation for the stream function and an advection equation for the concentration field. In fact the concentration equationfieldwill determine the position of the interface between thefluids.In this thesis, we will extend this model to its transient version. In this way, we will be able to study interfacial instabilities when a longfingerdevelops along the annulus. 8 Chapter 1. Introduction A longer description of the model will be given in chapter 2. This model was further studied by Pelipenko & Frigaard in [57]. They showed that in fact the model is well-posed. Also, for , certain fluid properties, a steady state solution exists, i.e. an interface that advances steadily along the annulus at the mean pumping speed. These solutions can befoundanalytically for concentric annuli and annuli with small eccentricities via a perturbation method. Their approach is also numerical, see [58]. They solve the system of equations using an augmented Lagrangian Algorithm, see [26, 30, 31]. This method fully describes unyielded regions of the flow if any. In [59], by the means of a lubrication approximation, Pelipenko & Frigaard deduce conditions, based on certain combinations of rheological and physical parameters,forwhen a good displacement takes place or not. They start from the assumption that a long finger has already developed along the annulus and assume pseudo-parallel flow away from the displacement front. By doing this they are able to predict the interface speed on both the wide and narrow side of the annulus and thus determine when the wide side interface moves faster or slower than that on the narrow side. This criterion defines whether or not a good displacement takes place, or whether the displacement front fingers along the well. 1.3.3 Displacements in a Hele-Shaw cell It is well known that the motion of a Newtonian fluid between twofixedparallel plates which are sufficiently close together (Hele-Shaw cell) is u = -^-^-zld2pdx y z) and v = -^-^-z(d - z) 1 2pdy y ' K (1.9) ' where z is the coordinate perpendicular to the plates and d is the thickness of the gap between the plates, see [8]. Averaging across the gap we get an expression for the gap-averaged velocity, given by: * = ~ v p (1-10) which is in fact the velocity field of a Newtonian fluid flowing through a porous medium of permeability k = yj. Thus, a direct analogy between the displacement flows in a Hele-Shaw cell and in a porous medium can be demonstrated. Note that (1-10) is Darcy's lawfora Newtonian fluid, we are interested in non-Newtonian flows, thus a nonlinear version of Darcy's Chapter 1. Introduction law has to be derived. In this section we present various models used in the literature. Hele-Shaw displacements for Newtonian Fluids We commence by describing the classical study of Saffman & Taylor, [64], regarding these type offlows.They studied the displacement of a viscous fluid by another in a Hele-Shaw cell. They found that the stability of the interface depends on the direction of the motion of the fluids, i.e. if the less viscous fluid is displacing the more viscous one, the thinnerfluidpenetrates the thickerfluidin the form of a long fingers, regularly spaced along the planar front. These fingers' themselves may destabilize and form sub-fingers, and the whole process may repeat leading to fractal patterns offingers.The displacement is ineffective and the interface is said to be unstable. This phenomenon has been termed "viscousfingering".If the more viscous fluid displaces the less viscous one, the interface will be stable. They also found that this is true whatever the density of thefluidsis, provided that the velocity is large. Since then, several studies can be found in the literature concerning Hele-Shaw displacements and we review only some of them, for an extensive review of the subject up to the mid 80's we refer the reader to [40]. Although thisfieldof work involves both immsicible and miscible displacements, we concentrate in the latter ones which are the type offlowsconsidered in the thesis. In 1986 Tan & Homsy, [74] studied the linear stability of a miscible displacement in a porous medium for rectilinearflow.They assume that thefluidshave the same densities and that the viscosity varies exponentially across the interface depending on the concentration of the fluid on top. Quasi-steady state approximation is considered because the dispersive concentration profile is time dependent. They showed that the linear displacement process is unstable for an unfavorable mobility ratio, i.e. viscosity of the displacingfluidis less than the viscosity of the displaced one. Thus the results are in agreement with the theory for immiscible fluids. Later, they performed numerical simulations, [75], to investigate the nonlinear effects of these instabilities. They found that once afingerbecomes large enough the tip of thefingerbecomes unstable and splits. Lajeunesse et.al., [46, 47], considered the downward displacement of a fluid by a lighter and less viscous one on a vertical Hele-Shaw cell. It seems that for a mobility ration 10 Chapter 1. Introduction greater than a critical value of 1.5, there exists a critical velocity for which the interface of the finger becomes unstable leading to a 3-D pattern involving regularly spaced fingers. Density driven unstable flows on Hele-Shaw cells have been also considered. Fernandez et.al., [25], studied density driven instabilities between miscible fluids in an horizontal Hele-Shaw cell, the fluids considered had equal viscosities and clearly, the heavier fluid was on top of the lighter fluid. They found that when the Rayleigh number, Ra, is small, the Rayleigh number being the ratio of convective and diffusive transport in the concentration (note that in the presence of net flow this number is just the Peclet number), diffusive effects are important and fingers with a wavelength i ? a _ 1 may appear. For the opposite case, Ra large, diffusive effects are negligible and a finger with a wavelength of the order of the gap width develops. Later, Goyal et.al., [36] conducted similar experiments with the difference that the fluids had different viscosities. They found that the instability is reduced if the viscosity is increased in either fluid. Another interesting finding was that the instability also depends on the thickness of the interface (diffuse layer). Hele-Shaw displacements for non-Newtonian fluids As stated before, in order to study these type of flows, we require a nonlinear version of Darcy's law, due to the fact that now the viscosity of the fluid is shear rate dependent. Barenblatt et.al. derived this relation for porous media flows, see [6]. They described the flow of a visco-plastic fluid through a porous medium by the following nonlinear Darcy's law with a limiting pressure gradient, s = -K -^)^ « = 1 0 |vp|>G |Vp| < G ' (i - n) (1.12) where G is the limiting pressure gradient defined as, G = CTyk' ' 1 2 with C being a dimensionless constant. In this 2-D system can be reduced to a single equation for the stream function if desired. Similar approaches to this type of nonlinear Darcy's law 11 Chapter 1. Introduction c a n be f o u n d i n the literature, see for e x a m p l e [4, 16, 34, 44, 53]. A c t u a l l y , the m o d e l d e r i v e d i n [10], w h i c h is the m o d e l we follow a n d e x t e n d i n this thesis, uses a n o n l i n e a r D a r c y ' s l a w w i t h a l i m i t i n g pressure gradient as well. H a v i n g stated this, we review some of the w o r k t h a t has been done c o n c e r n i n g n o n - N e w t o n i a n fluid displacements i n a H e l e - S h a w cell or a porous medium. W i l s o n , [84], considered the Saffman-Taylor p r o b l e m for a power l a w fluid d i s p l a c e d b y air. B y means of linear s t a b i l i t y analysis he found t h a t the g r o w t h rate of the i n s t a b i l i t y is 0 ( n - 1 / ), 2 where n is the power l a w index. L a t e r A z a i e z & S i n g h , [4] considered m i s c i b l e displacements of shear t h i n n i n g fluids by a N e w t o n i a n fluid or vice versa. T h e i r a p p r o a c h is the same as the one r e p o r t e d b y T a n & H o m s y , [74] for N e w t o n i a n fluids, the quasi-steady-state approach. T h e y found t h a t for the n N - N case, ( m e a n i n g n o n - N e w t o n i a n - N e w t o n i a n d i s p l a c e m e n t ) , the displacement c a n be unstable even for favorable m o b i l i t y r a t i o , this was expected due to the shear t h i n n i n g b e h a v i o u r of the fluid. F o r the N - n N case they found s i m i l a r results as for the N - N case. S t u d i e s i n w h i c h b o t h fluids are n o n - N e w t o n i a n are presented below. P a s c a l s t u d i e d the classical p l a n a r interfacial instabilities i n porous m e d i a for i m m i s c i b l e viscop l a s t i c fluids, [53, 54, 55]. H e used the M u s k a t approach, i n w h i c h one assumes a l o n g finger t h a t extends ahead of the m o v i n g interface, i f the v e l o c i t y o f the fluid w i t h i n the finger is faster t h a n the v e l o c i t y of the fluid outside, the interface is unstable. B y d o i n g this he showed t h a t i n contrast w i t h w h a t the theory of displacement flows i n a porous m e d i u m for N e w t o n i a n fluids says, even w i t h the presence of a n unfavorable m o b i l i t y ratio, i.e. less viscous fluid d i s p l a c i n g more viscous fluid, a positive difference of y i e l d stresses w i l l a i d the displacement. B y positive we m e a n t h a t the d i s p l a c i n g fluid s h o u l d have a higher y i e l d stress t h a n the d i s p l a c e d one for this to h a p p e n . H e also considered the effects of i n c l i n a t i o n of the m e d i a , a n d showed t h a t g r a v i t y forces t e n d to stabilize the interface as l o n g as the heavier fluid is under the lighter one. I n 1999 C o u s s o t , [16], considered the Saffman-Taylor p r o b l e m i n v o l v i n g two y i e l d stress fluids. H e found the same characteristics of instabilities as for N e w t o n i a n fluids except t h a t the w a v e l e n g t h of m a x i m u m g r o w t h c a n be s m a l l even at v a n i s h i n g velocities, i.e. i n s t a b i l i t y m a y o c c u r even i n the absence of g r a v i t y effects, as soon as the y i e l d stress of the d i s p l a c e d fluid is 12 Chapter 1. Introduction sufficiently large. Linder et al., [48], validated these results with their experimental work. 1.3.4 Interfacial instabilities There have been a number of studies of stability of immiscible iso-density multi-layer Poiseuille flows, in the two classical geometries of a plane channel and pipe. Whereas we are interested in the multi-layer Hele-Shaw cell , the interfacial instability is similar in both flows, so we will review of some of the most relevant work done for this subject. Both analytical and numerical studies of these instabilities have taken place. Explanations of physical mechanisms that govern this type of instability have been given in [38]. The earliest theoretical investigations of instability due to viscosity stratification is probably the classical study of Yih [88]. He studied the stability of two superposed fluids of different viscosities in plane Couette and Poiseuille flows, using an asymptotic method based on long wavelengths. He found that for equal densities of the fluids the flow can be unstable at any Reynolds number, however small. Hooper and Boyd, [41], produced the first studies of such flows in the short wavelength limit, for a plane Couette flow in which each fluid occupies a halfplane. Yiantsios and Higgins [87] consider a plane Poiseuille flow, in which 2 fluidsflowtogether down a channel. They extend Yih's results [88] for small wave numbers (long wavelength) to large wave numbers, accounting for differences in density and fluid layer thickness ratios, as well as the effects of interfacial tension and gravity. They show that the flow is unstable for large wave-numbers and is linearly unstable to a shear mode instability. In 2003 Ern et.al., [23], reported that instability exists in a shearflowof two superposed misciblefluidswith the presence of a continuous but steep variation of viscosity. They showed that Yih's results are valid even in the case of continuous changes in viscosity. Instabilities arise provided that the thickness of the interface is small and the Peclet number to be sufficiently large, i.e. close to the immisicible limit. Govindarajan, [35]foundsimilar results using a different profile for the viscosity variation. For the case of multi-layerflowsin porous media, Raghavan & Marsden, [62] showed, in the presence of density difference, if the heavierfluidwas in top of the lighter one, theflowwas 13 Chapter 1. Introduction unstable for all viscosity ratios, if the case is reversed, the lighter on top of the heavier, they claim that theflowwill be unstable only for large viscosity ratios (the more viscous on top of the less viscous). The isodensity case resulted stable for all wavelengths. Similar results were found by Gondret et al. in [33]. Although their problem is not the same as the one we consider in chapter 4, their results are fairly similar with the ones that we found in §4.5.1. Zeybek & Yortsos, [91], considered the parallel flow of two immisciblefluidsin a Hele-Shaw cell with equal densities. Through a weakly nonlinear analysis they found that the frequencies of the disturbances of the flow are strictly real and when the mobility ratio is different from one this theory predicts the presence of dispersive waves. They showed analytically and later experimentally the presence of solitary waves at the interface of thefluids.After the waves or solitons leave the channel the interface remained flat. These results reaffirm that without buoyant forces, the interface will be linearly stable. Gondret et al., [32] first considered the parallelflowof a gas on top of a viscousfluidin a horizontal Hele-Shaw cell. They presented a Kelvin-Helmholtz-Darcy theory which consists of Euler equations for inviscid flow coupled to Darcy's law for the viscous flow. They found that linear stability is governed by inertia only. Later, Hinch & Plouraboue, [60, 39], compared Gondret's analysis by considering the full Navier-Stokes equations and got similar results. 1.4 Overview As the brief account above indicates, unsteady displacements along an eccentric annular geometry often occur. Under this situation a long finger develops due to the fact that the displacement front will advance much faster on the wide side of the annulus than on the narrow side. When this happens, the interface between thefluidsis likely to become unstable. Despite the large number of theoretical, numerical and experimental studies, as shown in the previous section, the stability of multi-layer parallelflowsin a vertical or near vertical eccentric Hele-Shaw cell, to the knowledge of the author, has not been studied. The work presented in this thesis focuses on extending to a transient version the Hele-Shaw model for displacementflowspresented in [10] and to study the the stability problem just mentioned above. In addition, we investigate 14 Chapter 1. Introduction how a mud channel might be removed by pulsation of the flow rate. The thesis is in 4 main parts, which we outline below. 1. Extension of an existing Hele-Shaw type model to the transient regime, i.e. in which the velocity field satisfies a time dependent evolution problem. Fluid accelerations are considered, but not the nonlinear inertial terms. This is presented in chapter 2. 2. In chapter 3 we address theoretical questions for our model, i.e. existence and uniqueness of a solution, and various results on the continuity of the solution with respect to changes in the model parameters. Thus, we show that our model is well-posed. We also prove that the transient model decays to the solution of a steady state problem, provided that the concentration doesn't change. 3. Stability of multi-layer annular flows is addressed in chapter 4. Unsteady displacement fronts tend to elongate along the well on the wide side. Eventually, they become effectively multi-layer flows. We investigate if such flows are in fact stable. A wide range of analytical and numerical results on this problem are presented in chapter 4. 4. One particular multi layer flow is when the mud is not removed by the displacement and remains stuck in channel along the narrow side of the annulus. We investigate removal of this layer by techniques such as pulsation of the flow. This work involves mostly numerical simulation. The results are to be found in chapter 5. 15 Chapter 2 Model derivation In this chapter we give a physical derivation of transient model and justification for its usage. A brief outline of the chapter is as follows. In §2.1 we outline the pseudo-steady model of [10] for the bulk displacementflowand explain a number of phenomena for which that model is inadequate. This motivates derivation of the fully transient displacement model. In §2.2 we derive a reduced shearflowmodel from the Navier-Stokes equations and discuss a number of the modelling assumptions. In §2.3 we approximate the reduced shearflowmodel from §2.2 by a gap-averaged model. We show that this approximation is valid when the concentrationfieldchanges slowly. In §2.4 wefinallyderive the transient Hele-Shaw model from the gap-averaged reduced shearflowmodel in §2.3. This is the basic model that we study for the remainder of the thesis. In §2.5 we introduce a variant of the model in §2.4, in which the fluids are separated by an interface, as opposed to being represented by a concentrationfield.This variant is more convenient for studying interfacial stability. Since the concentrationfieldis only advected, (as is the interface), the two models are very similar. Their equivalence is discussed. 16 Chapter 2. Model derivation 2.1 The pseudo-steady model and some open questions In primary cementing the annular gap is comparatively narrow, compared to the other geometric length-scales. This is a large-scale, imprecise industrial process in which resolution on the scale of the annular gap is not always of interest. Therefore, the modelling approach that has been taken i n the past is to average across the annular gap, leading to a two-dimensional Hele-Shaw type model. This basic approach to modelling the displacement is described in detail i n [10] and analysed further in [57, 58, 59]. The pseudo-steady model of [10] relates to an entire annular wellbore, along which the geometry, may vary slowly. Here we will consider instead a single section of the well on which a constant geometry may be assumed. This is also the simplification considered in [57, 58, 59]. Assume therefore that fluid 1 (cement) displaces fluid 2 (drilling mud), along a uniform inclined annulus that is initially full of fluid 2. This situation has been modelled by the following elliptic problem: V-[S + f] a X(|V* |) TV 8 W s |V* |=0 5 = 0, (2.1) <=> \S \>jj, (2.2) <=* \S \<^, (2.3) s S where \& is the stream function. The unwrapped narrow annular space is (0,0 € (0,1) x (0, Z), s the annular gap half-width is H{<j>) = 1 + ecos7T0; e G [0,1) is the eccentricity, defined as A r / ( r — ri), see Figure 2.1 . The vectorfield / and the rheological properties of the fluids may 0 depend on time and space via the mean concentration of fluid 1, c(0, £,£). The function x is a positive increasing function of | V \ I / | , arising from the viscous shear-thinning behaviour of the S fluids, and Ty is the fluid yield stress. The concentration is simply advected: l[Hc} + -^[Hv c] s - W = H w + ^- [HW c} = 0, , = -Hv . z s s -j; 17 (2.4) a (2.5) Chapter 2. Model derivation F i g u r e 2.1: G e o m e t r y of the well T h e d e r i v a t i o n of (2.1), at fixed t, follows s t r a i g h t f o r w a r d l y from the reduced shear flow m o d e l : 0 = _ where r = 5 (T ^y,T £y) s S - dp ^ + d — T dp a m d m + g 6 , (2.6) d depends o n the v e l o c i t y gradients (2.7) (v ,w y) Sty S: at l e a d i n g order, a n d (<?<£> 5£) represents the g r a v i t a t i o n a l acceleration. These equations, plus the i n c o m p r e s s i b i l i t y c o n d i t i o n , are essentially averaged w i t h respect to y, that is, the r a d i a l c o o r d i n a t e scaled relative to the distance from the centreline of the annulus, a n d cross-differentiated,-to give (2.1). T h e pseudo-steady m o d e l above has been s t u d i e d i n some d e t a i l . I n a u n i f o r m i n c l i n e d eccentric annulus, i t appears t h a t the displacement is characterised by one of 3 different b e h a v i o u r s . a) If the interface between fluids advances steadily a l o n g the well at the same speed a l l a r o u n d the annulus, this leads to a g o o d displacement. T h e s e steady states have been observed i n the s i m u l a t i o n s of [10]. F o r s m a l l eccentricities, a n a l y t i c a l expressions c a n be found t h a t describe their shape, see [57]. T h e i r s t a b i l i t y has been characterised i n [58]. b) If a steady state displacement is not found, t h e n i n v a r i a b l y the interface r e m a i n s unsteady a n d elongates i n t o a progressively l o n g finger. 18 Chapter 2. Model derivation c) An extreme version of the unsteady displacement occurs when the fluids on the narrow side do not move, which represents a static mud channel. This can also be observed in the simulations of [10]. The lubrication model in [59] can be used to predict which of the above states occurs for any given set of process parameters and its fair to say that the mechanisms by which the finger and mud channel form are quite well understood. 2.1.1 Open problems and limitations The pseudo-steady model has a number of physical limitations and there are also aspects of the primary cementing process that are not accounted for in this model. 1. Since the velocity field is governed by a steady Hele-Shaw model for which the inertial and acceleration terms in the Navier-Stokes equations have been neglected, it seems unlikely that any type of interfacial instability can manifest, apart from the type of "viscous fingering" that results from gradients in mobility ratio. The unsteady displacement finger that advances on the wide side of the annulus can be interpreted as this type of instability. As the finger grows progressively longer, we find essentially a multi layer flow along the length of the annulus. Since these multi-layer flows are stationary solutions of the full pseudo-steady model, it is unlikely that there can be any fingering type of frontal instability of these interfaces. In reality however, it is well known that many viscous-viscous interfaces in shear flows suffer from interfacial instabilities, see e.g. [42]. In the context of annular displacements, Tehrani et al. report observing interfacial instabilities in their experiments as the interface becomes progressively parallel to the annulus axis, [78]. We study the onset of interfacial instabilities in chapter 4. When the interface is unstable this is likely to lead to azirhuthal mixing of the fluids. In the context of primary cementing, this effect could conceivably be either good or bad for the process. If the displacing fluid considered is a cement slurry, then the risk inherent mixing is that the slurry becomes contaminated and does not later harden when the displacement stops. On the other hand, 19 Chapter 2. Model derivation if the displacing fluid is a spacer fluid, then azimuthal mixing may lead to a better actual displacement of the mud than if the two fluids had remained separated by a distinct interface. 2. Although static mud channel formation is well understood, methods for removal of the mud channel are not well understood. One result of our study in chapter 4 will be to show that flows with an established static mud channel are linearly stable. Thus, interfacial instabilities do not appear to be a viable method for removal. We study this further in chapter §5. 3. Pulsation techniques have been advocated at different times for primary cementing, e.g. [19], with different claimed benefits. The effect of pulsation on displacement flows is not at all understood. To pulse the flow rate is perfectly possible with the pseudo-steady model (2.1), via the boundary conditions. However, we shall see in chapter 5 that there is little effect with the pseudo-steady model. The effects are quite different when a fully transient model is considered. 4. The effects of casing rotation and reciprocation on displacement are unknown. These are common practices in horizontal well cementing. We do not however study this problem, as it requires significant changes to the type of model that we present. Problems (l)-(3) require a transient model for the velocity field, since different timescales occur in the displacement. If the displacementflowsare still relatively slow, so that Hele-Shaw type averaging retains validity, then retaining the acceleration terms in the Navier-Stokes equations leads to the simplest transient models. 2.2 Modelling transient bulk flow cementing displacements As we can see in Figure 2.1, the geometry of the well can be modelled by an eccentric annulus, thus a cylindrical coordinate system (f ,#,£) is used to describe the well. From now on we will focus on only a uniform annulus. This will suffice to model a section of the well, for example, a 10m long casing stand or a laboratory experiment. Thus, £ is the measured depth, but measured 20 Chapter 2. Model derivation upwards from the bottom of the annulus. The inclination angle of the annulus is denoted by fi. The inner and outer radius of the eccentric annulus are denoted by, fj and r , respectively. 0 fV is just the radius of the inner cylinder and f is the radius of the outer cylinder, see Figure 0 2.1. We define the mean radius as f* = \(f + fj) and the mean half-gap width d = \{r — fj). 0 0 To model annular cementing flows we adopt the same approach as in [10]. We will concentrate on the displacement of drilling mud by a cement slurry, thus we have only one concentration field. In outline, we start from the Navier-Stokes equations, written locally on each section of the annulus. We observe that the annular gap scale is typically much smaller than either the circumferential or axial length-scales. This leads us to differential scaling of velocities and lengths, as is characteristic of Hele-Shaw flow modelling. The Navier-Stokes equations and concentration equation for the model are: „ „(du \dt A J Jdv . \ P ^ + u - V v ) Jdw 19... . 1d rdr r 06 1a,. ., = - - „ „ \ 0 { r \ 2 e f ) ^ rdr = d£ --r rd9 „ rgg n 6 6 + - , M - - dp or r 0. 19. ' y d 19. + lfl... , „ f ldp _ - ^ p ^ dp _ , n „. .... (2.9) • du, d£ where the gravitational vector is given by g? = gsm(3{€)cos6, | + + gg = g sin f3(£) sin 9, + where c € [0,1]. 21 g^ = gcosf3(£). = V.(6(c«)Vc) (2.12) Chapter 2. Model derivation 2.2.1 Non-Dimensionalisation The main aim of this section is to reduce the system (2.8)-(2.10) and (2.12) to the transient version of (2.6)-(2.7). We will average across the gap thickness using almost the same scaling arguments as in [10]. This is motivated by the fact that in a typical oil well we have the following situation, annular gap <C annular circumference <g; annular length-scale. Our goal is to derive a two-dimensional evolution model of the fluid motion, in azimuthal and axial directions. In contrast to [10] we will scale time with a viscous time scale and thus retain time derivatives terms at leading order but not the inertial terms. In the following sections we present geometry and velocity scales and the reduced systems of equations. Geometry and velocity scales Here we consider a single section of the well, of uniform geometry and length Z. We scale axial and azimuthal coordinates as: —, z= irr* <P = - 7T A measure of the narrowness of the annulus, 5*, is given by: 8* = ± f* The dimensionless eccentricity is given by: 22 Chapter 2. Model derivation Therefore the radial coordinate is scaled relative to the distance from the centreline of the annulus, The axial and azimuthal velocities are scaled with: U* = — A*' where Q* is the maximum flow rate and A* is the cross-sectional area of the annulus. The radial velocity is scaled with S*U*. We assume that the fluids can be described by the Herschel-Bulkley model, characterized by density, yield stress, consistency and power-law index. A characteristic scale for the rate of strain is: ' r*5*ir using this we can define scales for the shear stress, viscosity and pressure: f* = max[f ,y + k (Y) ], p.* = Z-, nk k fc k j* The dimensionless rheological parameters are defined by: and the fluid densities are scaled with 23 p* = ^f. 0 Chapter 2. Model derivation p*=max[p ]. k k Up to here the scales used to reduce our system of equations are the same as in [10], the only difference between our scaling and theirs is the timescale we use: _ p*{5*vr* f a p.* clearly this is a viscous timescale. 2.2.2 Reduced shear flow equations After the non-dimensionalisation we have the following reduced system of equations, at leading order in 5*: 0 = - 9 dy (2-13) 1 dp dv m d , = -7 d4> d- ^ p + T a dw + y dp m = -i p / d dy ge d-y ^ + T r a sT*' d<p s%' + . ( 2 1 4 ) . (2A5) <9£ where St* = ^ p*g* [r^*] 2 is the Stokes number which is usually <C 1. The boundary conditions at y = 0 and y = H are: at y = H, u = v = w = 0, u = T f y = T^y ( > =0, 24 at y = 0, (2.17) (2.18) Chapter 2. Model derivation The leading order rate of strain is characterized by: (2.19) 7 In regions where the fluids are yielded we have: 4m T V(i) r ~ (j) + 0(5*M iy V + 0(5*/V) (2.20) where 77 is an effective viscosity depending on the rate of strain and rheological parameters. We return to this reduced system of equations in §2.3, but first continue with the gap-averaged velocity and concentration fields. 2.2.3 2D averaged velocity field The aim of our modelling exercise is to reduce to a 2D model system of equations, from the shear flow model (2.13)-(2.16). To eliminate u from (2.16), we average across the gap and use the no-slip conditions at the walls of the annulus. A simplifying assumption made is that the annulus is assumed sufficiently narrow that the velocity components are approximated to 0(5*) by a velocity field that can be assumed symmetric about the annular gap centreline, y = 0, i.e. curvature effects only manifest at 0(5*). This allows us to average over only half the annular gap. The result of the averaging is: (2.21) where wdy. vdy Equation (2.21) is satisfied using a stream function formulation: 25 Chapter 2. Model derivation 2.2.4 Fluid concentration field The concentration of the displaced fluid is governed by (2.12). As our aim is to derive a 2D model for the large-scale motions of the fluids, we will want to work with quantities that are averaged across the annular gap. Thus, we proceed as in [10], by assuming that the concentration is uniform across the annular gap. This allows us to interchange the average of cu with the product of the the averages of c and u, in the convective terms. Such approximations are commonly made in hydraulic models, e.g. for the inertial terms, often including a prefactor to represent an assumed degree of correlation between the two averaged quantities. Here we have no a priori knowledge of whether the concentration and velocity fields are correlated and have taken the prefactor to be 1. With this assumption, we average across the gap, use the non-slip conditions at the walls and assume no flux of any fluid through the walls. The concentration satisfies: d (2.23) where the diffusivity has been scaled with a global diffusivity scale, D*, and the Peclet number is defined as: „ Trf*w* Pe = D* —. Typically, D* ~ 10" m s"\ thus Pe » 1 and typically Pe ~ 10 -10 . Denoting the scaled 10 2 10 12 annular length by Z, we may expect a diffusive layer of intermediate concentrations, in the azimuthal and axial directions, of thickness ~ y/Z/Pe -C 1. Therefore we neglect molecular diffusion in our leading order model. Thus, the concentration field will satisfy: ~[Hc} + ~{Hvc] + ^[HWc] = 0, (2.24) We shall assume that during the displacement, at the bottom of the well the concentration will 26 Chapter 2. Model derivation always be zero and at the top one, i.e. cement at the bottom and mud at the top. As well, natural boundary conditions at 0 = 0, 1 will be assumed. Note the only difference between (2.24) and (2.4) results from the different timescale used. We shall denote the ratio t* /t* where t* v av av = jff by e, i.e. ._ p*5* irr* U* 2 a A* Formally e —> 0 as 5* — > 0. However, although t* = 0(S* ), the ratio t*/t* is often 0(1). This 2 av is because 5* is small but not zero and -M^ is relatively large. In this thesis we shall consider P 'a both e — > 0, for which the concentration is fixed in time as the velocity evolves on the viscous timescale, and also e ~ 1. Since much of the complexity of the model is in the evolution of the velocity field, it can be helpful to isolate this evolution. Theoretical properties of the velocity problem (in the limit e —» 0) are studied in chapter 3. Later, in chapters 4 and 5 we consider the distinguished limit of e > 0 and S* —> 0. For this model, both velocity and concentration fields are fully transient and coupled. This allows us to study the transient phenomena cited before in §2.1. Apart from viscous and advective timescales, there may be timescales associated with the growth of any temporal instability or with pulsation of the flow. These may excite a response in either the velocity field or the concentration field. This is another motivation for retaining the time derivatives in both evolution equations Uniform concentration assumption The assumption of a uniform concentration across the annular gap is non-trivial. One way this might arise is via Taylor dispersion. In [85, 92] the authors showed that in the limit ePeG —• 0, where e is the inverse aspect ratio of the annulus and PeG is the Peclet number, (based on the gap width), diffusion effects dominate and a uniform or near uniform concentration field across the gap is achieved. For our problem, the Peclet numbers based on the gap width are in the range Pec ~ 10 — 10 and the aspect ratio in a typical well can not justify ePec <S 1. 8 10 27 Chapter 2. Model derivation Therefore, we are certainly not in the Taylor dispersion regime. Indeed, from the perspective of having large Pec, we might expect to be in the almost immiscible regime. Here either the displacements are effective and the uniform concentration assumption is reasonable, or otherwise long fingers are likely to form, with thin miscible diffuse interfacial layers of thickness Pe ^ . 1 G 2 Although such fingers are observed in the lab setting, we believe this scenario to be unlikely in the well for a number of practical reasons. First, along the well one generally has centralisers fitted at regular intervals. These will disrupt and mix the flow across the gap. Second, even without centralisers, the existence of a long finger extending stably over 10's or 100's of meters is implausible. There are various instabilities that may occur. For example,fingering& frontal instabilities were studied in [46, 47] in the context of the downward miscible displacement of onefluidby a lighter and less viscous fluid, between parallel plates. If frontal instabilities don't occur and afingerdevelops, there have been many studies of interfacial instabilities in shear flows, see e.g. [5, 42, 67] for an overview. In general, suchflowsare not robustly stable to all wavelengths and rheology combinations. The situation is not apparently helped by having a thin diffuse interface layer over whichfluidproperties change. Recently in, [23] the authors studied miscible multilayerflowswhere the viscosity varies rapidly between two constant values. As the layer thickness decreased, they showed that interfacial instabilities arise in the same way as for the immiscible fluid limit (with no surface tension). Thirdly, for the cases that may be stable in an idealised mathematical or lab setting, we must recall that the annulus outer wall is perhaps rough with small scale unevenness, that the pump pulsations will naturally occur, and that there are also likely to be azimuthal secondary flows. It is almost impossible to judge the effects of all the above effects, and our assumption of an homogeneous concentration is really just a statement that our model is not designed to study small scale flow features on the gap scale. 28 Chapter 2. Model derivation 2.3 Reduced shear flow models At this point we would like to proceed in the same way as in [10] andfinda transient displacement Hele-Shaw type of model. Unfortunately, the model (2.13)-(2.16) does not lead to'a closure relation between the modified pressure gradient and the gap-averaged velocity. Thus, it is only a simpler 3D model of the flow. To circumvent this difficulty we develop an approximation to (2.14) & (2.15) in which the solution remains close to the solution of the transient model, i.e. the difference between the solutions of both models is O(e); and allows a closure relation between pressure gradients and the velocity. Our starting point for the model derivation is the unsteady leading order shear-flow model, (2.14) & (2.15), which we write as: put = (2.25) G+^-r, dy where u = (v, w) and G = (G G ) EE h + , - | + g^j . S (2.26) g$ The system (2.25) is supplemented with the conditions (2.13) & (2.16) and (psin/3sin7T0 pcosjA sf (9*>9t)={ ' (2 27) As well as the system (2.25), which we shall refer to as the transient system, we shall consider the steady system, (2.6) & (2.7), which we write as: 0 = G + —T . dy s (2.28) S The subscript s is used to indicate that we are seeking a velocityfieldu = (v ,w ) that s s s satisfies (2.28), where potentially G ^ G, (e.g. because thefluidconcentration will be advected S differently by u than by u). As with (2.25), the pressure and concentration are independent s of y. In addition to (2.28) the incompressibility condition (2.16) is satisfied: ~ dy + d<j> 29 + Chapter 2. Model derivation Finally, we consider the gap-averaged approximation to (2.25): . pu , = G +^ T , a t where again the subscript a a (2.29) A denotes that we seek a velocity field u = a (v ,w ), a a and that potentially G ^ G; (2.29) is supplemented with the incompressibility condition: a du a dv dw a a The overbar indicates a gap-averaged quantity, e.g. u = a V a = H J V a 0 y ' W a H j = (v ,w ): a W a a ' V It is important to note that the time evolution in (2.25), (2.28) and (2.29) is distinctly different. Starting with the steady problem (2.28), note that given G at a point ((f), £), we are able to find 3 a unique solution (v ,w ), (and could then integrate the incompressibility condition to find s s u ). a All the time dependency in (v ,w ) must come from the vectorfield G , which is determined s s s from (2.1), with appropriate boundary conditions, and the advection of the concentration field. When G and c vary only slowly with time, (v , w ) are (pseudo-)steady on the 0(1) timescale. s s s s Considering now (2.25), we see that this is simply the parabolic extension of the elliptic system (2.28). If G and c are time invariant, (v, w) will still evolve with time, i.e. (2.25) is an evolution problem for (v,w), regardless of G and c. Turning now to (2.29), we could interpret this as an integro-differential equation for (v ,w ), a a but it is slightly unclear if, for specified G and a c, a the evolution of (v ,w ) is properly determined. Rather, given time derivatives of the averaged a velocities (v ,w ,t), att a w a e c a n determine unique (v ,w ), (i.e. this is the same problem as (2.28) a a at fixed t). It remains to determine whether (v,w) evolve in a manner consistent with these solutions. Apart from retention of the time derivative terms, the model and assumptions leading to (2.25) are identical with those in [10]. Inertial terms have been neglected under the asymptotic assumption, e <SC 1. We shall be interested in temporal changes in the annulus that occur over the slow timescale e . It will however be the gap-averaged system (2.29) that forms the basis _1 30 Chapter 2. Model derivation for the generalisation of (2.1), that we consider later in chapter 3. Since (2.29) is not derived directly, we should like to estimate how close the solutions of (2.25) and (2.29) are. We consider first the case n = 1, in which case the fluids are Bingham fluids. Later we discuss the case n < 1. 2.3.1 Stability of u and u , at (d>, £): 8 For n = 1,'problems (2.25) and (2.28) are of a form considered by [20, 30]. Multiplying (2.28) by u. and integrating with respect to y, we have: s a (u , s s = L(G -u ) u ) + j (u ) s s s s (2.31) s where (u,v) s = K / Jo Jo Uy-Vy S (2.32) dy, f" j (v) = TY, / \v \ dy, Jo L{f) = [ fdy. s S (2.33) y (2.34) Jo From this we can deduce that: \us\\m<—\\G \\ ,, K TT s ' L (2.35) S from the Cauchy-Schwarz and Poincare inequalities, where ||M ||i/i is the semi-norm: s 1/2 Mm = ( Jo/ N 2 y d which is equivalent to the H -norm. Using the methods in [20, 30] we may assert that there 1 exists a unique solution u . Time-dependency enters u only through G and c . We assume s s s s that u is differentiable with respect to each of the physical fluid properties, which depend on s c , and with respect to the components of G . Therefore, we may suppose that s s u Stt = eb(et), for some vector b. 31 (2.36) Chapter 2. Model derivation Now we consider stability of (2.25). Multiplying (2.25) by (u — u) and integrating with respect s to y we have: p u - {u —u)dy t > s L(G.(u — u)) — a(u, u — u) s s Jo (2.37) -j( )+j(u), Us where a and j are defined analogously to a and j , but with non-subscripted rheological s g constants. Combining (2.36) with (2.28), we have: pu = Sft G + s (2.38) + epb. T, s y Multiplying by (u — u ) and integrating with respect to y gives the following inequality: s f H p / u Jo Sit • {u - u ) dy > 3 L(G • (u - u )) - a (u , S s s -js(u) + j (u ) s s s + epL(b • (u - s (2.39) u - u) u )). s Summing the two inequalities we have o"57 / z at J \ s~ \ &y u u 2 < -a{u - u,u s s - u) + L((G - G ) • {u - u )) s s 0 +epL(b • (u - it)) + ^1 - -^j a (u , s +(l-^) s u - s u) s [ ? - ( « ) ( 2 - 4 0 ) Using the Cauchy-Schwarz inequality we have: | ^ | | U s - « | | 2 a < - \\u,-u\\ K (2.41) 2 H i + (\\G-G \\ 2 + + (\K - K\ \\U \\ I s L S H S ep\\b\\ 2)\\u -u\\ 2 L s L + \TY, -TY\H ^ 1/2 S \\U S - u\\ i. H We may split the first term on the right hand side in two. From the Poincare inequality: \\u -u\\ i>^\\u -u\\ 2. s H L s Therefore either we have that: 2H (\K -K\ \TY,s-rY\\ S IK - < — \-^r\\ '\\» G 32 + gi/2 ) . / 0 4 0 x ( - ) 2 4 2 Chapter 2. Model derivation or we have the decay equation: d p— \\u -u\\ 2 s K7T - ^ \ \ u - u \ \ 2 < L s L 2 + (\\G-G \\ 2+ep\\b\\ 2), s L L from which, for t € [0, T] |« -u|| 2(t) s L \G-GsWv < / Jo -Xi(t-s) + ep||b|| 2 e d g L +||t -n|| (0)ets Alt L2 , (2.43) with Ar= • f ^ ) lin i n ^ ^ J . m 2 te[o,T] In either case, we note that locally 1\u — M|| jji is bounded by the difference in the data between s steady and transient problems. 2.3.2 Stability of u and u , at (tp, £): a a We proceed analogously to the above, with the variational form of (2.29) we modify, to include the term pHu , s t u t, s< • (u - its) > a L(G • (u - u )) - a (u , S s a • (u - u ) s a > a s +j (u ) a s s s s a (2.44) s + epHb • (u - u ). s a L(G • (u - u )) - a (u , a u - u) n a a s u - u) s a -ja(u )+ja(u ). s (2.45) a Summing these two inequalities, we have 2^^ll s-«o||ii M The latter i.e. -j (u ) pHu ,t & (2.28). < * -a (u -u ,u -u ) + a +L((G s S + a s a epHb-(u -u ) s a - G ) • (u - u )) a s ^l-^ja (u ,u s s a a -Us) + (l-^j[js(u )-js(u )}. a 33 s (2.46) Chapter 2. Model derivation Again we split the decay term in half, and bound above using the Cauchy-Schwarz and Poincare inequalities. Either we have: ^H\\G \\ 2 S < \\U -U \\ i s a H \K -K \ a a . h \TY,s ~ TY, \ L K a + \\G -G \\v> a 1 2 a 7TK K s .2H ' a (2.47) — , 7TK a a or we have the decay equation: p d ,. K TT . 2 a r from which, for t G [0, T] ||« -«a||Li(t) 8 < f eH\b\e- ^-^ Jo x As ' +IK-«alb(0)e- A 2 t , (2.48) with , f H 7T . A = mm i 2 t€[0,T] \ 2 A 0 }. \2H pj 2 Again we see that ||u — w||x,i remains uniformly bounded. Note that we cannot expect results s on \ \u — u \\ 2, since the averaged problem leads only to an equation for evolution of the mean s a L velocity. Summary and further comments Combining (2.43) & (2.48), we see that ||tt — 0 u|| i(t) L is bounded by differences in the data between the steady, transient and averaged problems. When these differences are uniformly 0(e) over the slow timescale, we can certainly expect that ||« — u|| i(t) = 0(e). Since we , a i work with the gap-averaged velocity fields, this is the type of approximation result that is needed, i.e. the averages are close provided e is small. In words, this means that provided the flow has sufficient time for viscosity to act, before the fluid properties or other process features change locally, then it is valid to interchange the time derivative and averaging operators. This is valid in the sense that the averages of the velocity fields, taken across the gap, remain close. More generally, we consider boundary conditions on the annulus that vary only on the slow timescale e , on which timescale c also evolves. Decay of initial errors occurs on an 0(1) -1 34 Chapter 2. Model derivation timescale, after which the errors are simply due to the differences in G and the material properties. These are governed by the evolution of the gap-averaged velocity fields, which takes place on the slow timescale. We should therefore expect \\u — u\\ i(t) = O(e) on a timescale a [0,T], where T~ ~ o(e), but T l - 1 L ~ e is unrealistic. The exact timescale will of course depend on the particular flow studied. For example, for a stable single fluid flow, one would expect \\u — a = O(e) on the slow timescale. If we wish to consider instability of a steady solution, we must expect that the model breaks down as soon as 0(1) differences in the data evolve, i.e. depending on the growth timescale for the instability. The arguments for n < 1 are more difficult, especially for two-dimensional flows. To circumvent this, we may note that the velocity u, governed by (2.25), decays to a vector in the direction of G on the 0(1) timescale, (i.e. the component orthogonal to G decays exponentially). Thus, we may consider a one-dimensional problem. For such a flow, we can proceed essentially as above. On summing the variational inequalities we use the bound: for C = C(n,H); (see Zeidler [90]). Thereafter, we may use Holder's inequality to bound \\u — v\\ i terms and the Poincare inequality on W (0,H) 1,n+1 H 2.4 for the rest. Eccentric annular Hele-Shaw displacement transient model Finally, we are in a position to derive the gap-averaged transient model. We commence with a characterization of the fluid rheologies. 2.4.1 Rheological assumptions The fluids are considered to be shear thinning generalised Newtonian fluids, and typically have a yield stress. After the Hele-Shaw scaling, the principal components of the rate of strain are i<t>y ~vv a n d 7£y ~ y w Thus, 7~[«2+«fl . 1/2 35 Chapter 2. Model derivation The leading order shear stresses and rates of strain are related by a law of form: nj=r){i)iij, '• ij = 4>y,t,y, ( -49)' 2 where 77 is referred to as the effective viscosity. When the fluid has a yield stress, 77(7) —> 00 as 7 —> 0, and T{j becomes indeterminate. In [10, 57, 58, 59] it is assumed that the fluids are Herschel-Bulkley fluids: 7Vi — 7 = 7ij. 1 T > TV, 0, (2.50) T < 7 V , where «, n and T V are the consistency, power-law index and yield stress of the fluid, respectively, (which will depend on c). Since the fluids are shear-thinning, we consider 0 < n < 1. Here we modify the Herschel-Bulkley model slightly, and will replace (2.50) with 7 . 7 = lij, T>Ty, (2.51) 0, r<7v, where p^ is referred to as the high-shear viscosity. There are 2 motivations for this modification. First from the physical perspective, use of (2.50) implies an effective viscosity: *?(7) = I T ' 1 + —, 7 and therefore 77(7) —> 0 as 7 —• 00. This is unrealistic, since suspension viscosities commonly approach a Newtonian plateau for large 7, i.e. the Herschel-Bulkley model is really intended to model low-shear behaviour. The second motivation is purely mathematical. Without the highshear viscosity we are later forced to work in subspaces of L 1+n and W , rather than I? and l,1+n H . Although most of what we prove in the Hilbert space setting can be extended to L 1 1+n and W ' , it is less convenient and involves additional analysis. From the numerical perspective l 1+n this is anyway unnecessary, since numerically we end up working in finite dimensional subspaces, which lie in both {L , W^ } and in {L , H }. 1+n 1+n 2 1 In practice we shall suppose that p.^ is small, by which we mean that at typical shear rates, 36 Chapter 2. Model derivation (which are Os(l) due to the scaling), we expect Moo < KJ In this sense we may regard p^ as a small regularisation parameter. We emphasise however that this is not the common regularisation that is used frequently to avoid proper treatment of the yield stress. Here it is vital that the yield stress behaviour is retained, since occurrence of staticfluidregions in the annulus is commonplace in cementing. Fluid concentration effects Finally, if we have intermediate fluid concentrations occurring in the annulus, it will be necessary to specify the physical and rheological properties of the mixture. We consider this purely as a problem of model closure. For simplicity we use linear interpolation of the physical properties. Whilst this is certainly not realistic for allfluidscombinations, there is no universally accepted closure. Nothing that follows in the model derivation or analysis depends explicitly on the closure law for the mixture, i.e. provided the closure is reasonably smooth and therefore measured mixture properties could be utilised later. The numerical simulation results evidently will depend on the mixture law adopted. TY(C) = TY I 2 C + 7-^1(1 - K(C) = c), n(c) = n c + ni(l — c) K2C + p(c) = P2C + 2 KI(1 - c), pi(l -c). (2.52) (2.53) 2.4.2 Hele-Shaw model derivation Wefirstintegrate (2.16) across the annular gap to eliminate u: which prompts definition of a stream function \I/(c/>, for the gap-averaged flow: (2.54) 37 Chapter 2. Model derivation Turning now to the system (2.29), we observe that the pressure and gravitational terms do not vary with y. Therefore, integrating twice with respect to y we have: Ua(y,<f>,£,t) y ^ v(y) I [G -pU,a,t] = a (2.55) dy, where u = (v, w) and dp G Therefore \G a u a — pu ,t\, a psmBsinircf) "d~4> a is instantaneously parallel to G a at each fixed dp St* + pcosB '~"di Writing - pu ,ta (2.56) St*~ and s(y) = \u \{y) a A = the speed s is related to the modified pressure gradient (<fi,£,t), A via the one-dimensional shear flow problem: ds dy dy A ds dy For generalised Newtonian fluids this problem is straightforwardly solved, either analytically, or numerically by simple quadrature. However, we are interested in the gap-averaged speed: = / s(y) dy = - / Jo Jo j r-H = yi{y) y J d Q r-AH ^2j = ^( ) T T d - (2.57) r where from (2.51), for r > ry, and ( T ) is obtained by inverting this monotone function. We may observe straightforwardly 7 that (r) is a strictly monotone C°° function of r > ry. If n < 1, then as T —> ry, (r) — >0 7 7 and the Herschel-Bulkley term dominates: — T 7 i 1/n Ty * 1 As r —> oo, the high-shear viscous term dominates: 7 T — TY Mo In between these limits numerical inversion is needed to define J(T). If n = 1, then T — 7 Ty : Poo & 38 Vr > ry. Chapter 2. Model derivation The expression (2.57) defines the closure relationship between the gap-averagedflowrate | V * | and the modified pressure gradient A. Evidently, as A —> ry/H the range of the integral vanishes and | V * | 0. Since the integrand, TJ(T), is strictly positive, increasing faster than linear and is C°°, for r > ry/H we see that |V\&|(.<4) is also C°° and increases strictly monotonically, for A > ry/H. Inverting, this relation, we may write A = A(\W\). For |V*| > 0 we have that AdV^I) is strictly positive, (bounded strictly below by ry/H), and strictly monotone. It is convenient to separate the yield stress effects. Below, we shall write A(|W|)=x(|V*|) + I £ . Thus, x(|V*|) represents the part of the modified pressure gradient surplus to that required to overcome the yield stress locally. We may rewrite (2.57) as: W " WT^jHfl <' > 2 58 Although we have focused on the relation between | V * | and A, (equivalently x), we note that via the constitutive law and the limits on the integration in (2.57) & (2.58), we have the following parametric dependency of these functions. x= | W | = IWKxjtf.TY.K.n./Zoo), x(\W\;H,Ty,K,n,p, ). 00 Returning to (2.55), we average across the half-gap, y £ [0, if], to give: j" -faiyty, *. = [G -pu ±[ a aA from which we see that u is also instantaneously parallel to the vector [G a <9# / <9# dz'd<t>) 2 dp d +94 V 2 — pu ,t]' a dp e<i> " d4>dt d$ \ d£dt p (2.59) +9i p A |W| For | V * | > 0, replacing A with x + ry/H, this implies that <9* dp 2 psmfSsmiTcj) d£dt d<p St* _dV dp pcos/3 _ d(f>dt d^ St* 2 p 39 _ (x(|V$|) + ry/H\ \ |V*| J di [x(\V*\) + TY/H\dV V J d(f> <9# (2.60) Chapter 2. Model derivation Cross-differentiating to eliminate the pressure, wefinallyarrive at: V • [pV*] - - v • t Yx(|v*|) + 7 Y / f f LV iv*| (2.61) v* + f /here / . cos/3 . , sin/5"sin7Tc/>\ f=(p(-c)^^)-V^js If | V*| ,„ „„x (2-62) = 0, we may still cross-differentiate to eliminate the pressure, except that the right-hand side of this system is indeterminate. We may write: V-[pV*t] = -V-[S + f], (2.63). where S= ( x ( | V ^ 7 Y / g | ) V* o |V*|=0 |S|>,y/2f, (2.64) \S\<TY/H. (2.65) We note also that: d4> depot 2 dp pcosB <9£ St" d^ dp ^<9£(% <9c/> 2 ' p sin B sin Si* Equation (2.63) is the classical formulation of the evolution problem for nq (2.66) (effectively giving the gap-averaged velocity in the annulus). We shall consider this in a more rigorous setting below in §3.1. For some intuition, observe that if n = 1, we consider a constant concentration, zero yield stress and assume that the annulus is concentric, then (2.63) is simply: pA^ t = -3(K + Moo)A#, with suitable boundary conditions and initial condition. 2.4.3 Boundary conditions Equation (2.63) is supplemented with the following boundary conditions: •*((U,t)=0. 40 (2.67) Chapter 2. Model derivation tf(U,t) = Q(*), (2-68) on the wide (<> / = 0) and narrow (0 = 1) sides of the annulus, respectively. Here Q(t) = Os(l) represents the total flow rate through the annulus, appropriately scaled. Conditions at the ends of the annulus are harder to specify, and depend largely on the situation that we are modelling. In general we shall suppose that large variations in the fluid concentration occur away from the ends, i.e. we are interested in displacement phenomena away from the ends. If we consider our constant geometry section to be a section of the well, then appropriate are: d$> S = O^^(0,O,i)=O. 5 • (2.69) a* 5 = 0 ^ — (<f,,Z,t)=0. 5 (2.70) Alternatively, if we model a lab-scale pilot experiment, we may impose e.g. the uniform inflow condition *(&0,*) = Q(i), (2.71) in place of (2.69), retaining the outflow condition (2.70). Finally, if we consider Z » 1 so that the flow close to the ends of the annulus is far from any concentration variations, we may calculate appropriate one-dimensional flows at the ends, which correspond to stream functions: *m(0) at z = 0, and * u t ( A ) at z = Z, respectively. We might then impose < 0 *(>,(),*) = * (0,t), (2.72)' tt^.Z.i) = *„*(>,*), (2.73) in in place of (2.69) k (2.70). The Dirichlet conditions (2.72) k (2.73) are easiest to handle analytically, and we assume this below unless otherwise stated. The system (2.63) with boundary conditions (2.67), (2.68), (2.72), (2.73), and equation (2.24)for the concentration will be the main object of study of the thesis. In the next chapter we will establish theoretical results that relate to the well-posedness of this model. In chapter 5 we solve this model numerically to study static mud channel removal. Chapter 4 addresses interfacial stability. For this study a fluid concentration is not as easy to work with as is an interface. Therefore, a variant of the model is derived below. 41 Chapter 2. Model derivation 2.5 Interface tracking model derivation Later in chapter 4 we consider a version of model (2.63)-(2.65) in which the fluid is separated cleanly by an interface into twofluidsof different constant physical properties. Below we derive this model. Derivation and jump conditions at the interface From §2.4.2 we have, V-[pW ] = - V •[£ + /] (2.74) t Let us assume for the moment that * e C (il). Now, let us multiply (2.74) by a test function 1 v such that, v :R 2 x [0, oo) -> R is smooth and with compact support on fl. After multiplying (2.74) by v and integrating by parts, we have 0 = / ( V * i + S + f)Jo. P Vvdfl = - / V • (pVtf* + S + f)vdfl (2.75) JQ We derived this equality assuming that * £ C (fi), but note that (2.75) has meaning even if 1 is only bounded. Suppose now that ^ is smooth on either side of a smooth curve C dividing fl into Qi and fl2, where £l\ is that part of fl on the left of the curve which has onlyfluid1 on it and 0 that part 2 on the right of the curve withfluid2 on it. Suppose that ^ is a solution of (2.74), as we will prove in chapter 3, and assuming that itsfirstderivatives are uniformly continuous in fl\ and 42 Chapter 2. Model derivation fl2 we can choose a test function v with compact support infii,thus (2.75) becomes, 0 = / (pWt + Si = - / V • (pVtf* + fi) • Vvdfl (2.76) + Si + fi)vdfl, Jili JUi where the subscript 1 means that we are considering the properties offluid1 only. Because fli is full offluid1 then V • fi = 0, thus PiA$ = - V • Si in t fii. (2.77) Likewise, p At> 2 t = -V• S (2.78) mfl 2 2 Now let us select a test function v with compact support in fl, but which does not necessarily vanish along the curve C. Then, using (2.75) we have 0= / (pV* + S +/) • t Jn Vvdfl (piVtft + Si + / 1 ) • Vudfii - f Jni = + [ (p V* + 2 Jih t S 2 + /) • 2 (2.79) \7vdfl . 2 Now, since v has compact support within fl, we have 0 = / (piW + Si + fi) • t Jn Vvdfl! = + = - / (piA* + V • (Si + fi))vdfl t Jni f {(piWt + Si + [ {(pi^t + Si + fi)- }vdl Jc x fJ-n^vdl. (2.80) ni Jc because of (2.77) and ni denotes the outward normal to C from fli to f2 . 2 Similarly we have 0 = / (p V* + S + f ) • Vudft 7n 2 t 2 2 2 = - / (p A* + V • (S + / ))i;dr2x Jfi + / [(p W + S + / ) • n }vdl 2 t 2 2 2 Jc = Adding (2.80) and (2.81) we have 43 2 t 2 2 2 / [(p V* + S + / ) • nzjvdi Jc 2 t 2 2 (2.81) Chapter 2. Model derivation (p Vtf + S + f ) • n\\ = 0, fc t k (2.82) k along C. Thus, our interface model can be written as follow, p A* + V-S fe t fc = 0, (^Oefifc, fc = l,2, (2.83) where the index k is used to denote the fluid. Equation (2.83) is supplemented with the following boundary conditions: tf(fU,t)=0. (2.84) *(l,^t)=Q(t), (2.85) on the wide (0 = 0) and narrow (0 = 1) sides of the annulus, respectively. Here Q(t) represents the total flow rate through the annulus, appropriately scaled. For much of this section we shall assume an infinitely long annulus, with suitable conditions on \& as £ —> ±oo. The two fluid domains, fifc,fc= 1,2, are separated by an interface, i.e. the smooth curve C, that we represent as the following level set: F(0,£,i)=O. (2.86) The function F(<p, £, t) satisfies a gap-averaged kinematic equation: HdF 7 a r d*dF r - ^ a 0 WdF + 7 3 0 5 ? = 0 ' ( 2 - 8 7 ) where H(<f>) is the half-gap width and e denotes the ratio of viscous to convective timescales. At the interface the following two jump/continuity conditions are satisfied, *| = 0, (2.88) VF| = 0. (2.89) 2 p cos/3 p sin/3sin71-0 fc , fe The first of these is simply continuity of c 2 whereas (2.89) gives the jump in normal derivative of*. 44 Chapter 2. Model derivation Equivalence with concentration model For simple and suitably smooth interface configurations, the kinematic condition for the interface can be related to the concentration equation, as follows. Let us assume that for each £ and t we can define <!> = {<l>M(Z,t):c(<l>M,Z,t) = 0.5} Assume that the region where the concentration field changes smoothly from 1 to zero is O(Sd), see Figure 2.2, where 5d is the diffusive layer thickness, which is of Oil/Pe / ). 1 2 We know from [85, 92] that if we consider the limit Pe —» oo we approach the immiscible limit of our model, and we can consider our domain to be separated cleanly by an interface into two domains which contain fluids of different constant physical properties, see Figure 2.2. Thus we can define our interface as 0.5 0 J± • 1 • Figure 2.2: Example of concentration field at fixed £ and t. Note that if we assume that the fluids are separated cleanly by for 4> G [0, <S>i{(,,t)) and V therefore 45 t), we have that c(<f>, £, t) = 1 Chapter 2. Model derivation Jo. and ail ^w=y -^r^m-^ 0 (- ) 2 at <{> = 91 0M, (2.92) Subtracting (2.90) and (2.92) we have / r<t>M \ e d0+ — / ^ ^ = - — ^ +— +0^ at^ = 0 . M 2.93) Similarly, we can get Integrating (2.24) with respect to <j> from 0 to C/>M, taking the limit Sd —> 0, i.e. the immiscible limit, and substituting (2.93), (2.94) and (2.95), we get H dfa <9# a* dfa -e ^ OT+ ^oi 7 7o<p ^ : ^ oi F + n = 0 ± at = , Therefore a natural choice for F(<j>, £, i) defined in (2.86) will be, F{fai,t) = j -fa{i,t). ( ) 46 ,„ 2.96 Chapter 3 Qualitative results Mathematical validation of our model is presented in this chapter. We prove that the evolution equation (2.63) is well-posed. For the Herschel-Bulkley model, regularised by the high-shear viscosity discussed in §2.4.1, we can work in a Hilbert space setting. This provides a considerable simplification to the analysis. Now, we present an outline of the chapter. i) In §3.1, first we present some preliminary results on the behavior of x d ^ * ! ) a n d * which will be needed for the existence and uniqueness proof. Secondly, we define the variational problem associated with the system (2.63), we also define the subgradient of the corresponding functional for the minimization problem (3.26). Finally, we follow closely the results found in [24] (section 9.6). The section referenced introduces the theory of certain nonlinear semigroups, generated by convex functions; theory also known as gradient flows. Using these results we prove that a solution for (2.63) uniquely exists and that: *GC([0,oo);fr (fi)) 1 > * €L ([0,co);ir ("))oo 1 t ii) In §3.2 we consider theflowof afluidwith a fixed set of physical and rheological properties along the annulus, i.e. these may depend on (0, £) and a fixed concentration field, i.e. e — + 0. We assume that one of the fluid properties undergoes a continuous change everywhere in the annulus, for example from a density field p\ to a densityfieldP2, while the other properties remain unchanged. We present various results on continuity of our solution with respect to the rheological and physical properties of the problem. The results for 47 Chapter 3. Qualitative results the transient model are summarized below, 1. Let p\ —> P2, thenfort G [0,T] | | p i - P2||i«»T + | | V ( * i - * )||L*(0) > ^ | | V ( * i 2 tt )||£ 2 2 9 (3.1) 2. If 7 y i —> 7 y 2 , we have liny - r y | | 2 T + ||V(*i 2 L *2)||L*(0) > C ||V(*i 2 * )HL2 2 (3-2) 3. Finally if K\ —»/c we get 2 I K - « 2 l U a r + | | v ( * i - * ) H L * ( O ) > c ||v(*i - * 2 )||| 2 2 2 (3.3) where the constants do not depend on the solutions or the rheological parameter that varies, i.e. we have continuity with respect to the initial data and physical parameters. iii) In §3.3 we derive a number of qualitative resultsforthe steady and unsteady problems. These results do have a physical interpretation, which we give, but also serve as a useful test in verifying that our numerical procedure is correct. We show that the total energy dissipation rate in the annulus increases with any of the yield stress, consistency or high shear viscosity. Also we present the decay rate of the solution for the transient model to the solution of the steady model, taking into account that we have a fixed concentration field, i.e. e —* 0, ||V(* - * )|| 2 S L < ||V(* - *s)\\i?Q)e- Kt where K = 3.1 1C Pmin for C<1. Existence & uniqueness of * Here we shall prove existence of a solution *(0,£, t). We are concerned primarily with the evolution of <I> in t. According to (2.24) the concentration c evolves slowly on the timescale e . _1 Consequently we will regard c as fixed in time, (i.e. formally we consider the limit e —• 0), saving 48 Chapter 3. Qualitative results for future work the analysis of the fully coupled system. Thus, our focus is on establishing that we have a well defined evolution problem for The annular flow domain is fi = (0,1) x (0, Z). We state first some physically motivated assumptions, that we will assume throughout. A l The concentration c(</>, £) £ L°°(fl), and is bounded by 0 and 1. The physical properties of the fluid are all smooth functions of c; p^, K, p and n are strictly positive, ry is semi-positive. All are bounded above and in particular n < 1, (shear thinning fluids). A2 Theflowrate through the annulus, Q(t), and the various pressure gradients in the annulus are bounded. With Al above, this implies that (S + pV*t) e [L°°(ft)} , 2 since we have that , S + pV* 0 t = ( , dp pcos/3 dp psin/Jsmmft _ _ _ ^ _ , _ _ ) (3.4) A3 The annulus eccentricity e satisfies 0 < e < 1, which implies that 1 + e > H(4>) > 1 — e > 0. Since we often consider a displacement between 2 fluids, with intermediate concentrations confined to the interior, at times we shall also assume the following. A4 As £ — > 0 and £ — > Z, the physical properties of thefluid,TY, p^, K, p and n all approach constant values. 3.1.1 Preliminary results We commence with a number of preliminary results. Properties of %(|V*|) P r o p o s i t i o n 1 For monotone; |V*| > x(|V*|) -> 0 as 0 we have that |V*| -» x(|V*|) is C°°, strictly positive and strictly 0. Proof This follows from the properties of A(|V*|). See the discussion in 49 §2.4.2. Chapter 3. Qualitative results Proposition 2 The function x ( | V * | ) is bounded below by XN(\^^\), defined implicitly for |V*| |V\IJ| > 0 as follows: (3-5)' 3Moo g X /A j B ( X B + 3 * andxHB{\^^\), = Z?*L, |V*| = V XB(|V\I>|) = J 1.57y/ff) ;3/ioo . ^( X' B/ +_ Tr y// F, )/ . (3-6) 2 7—"TON 7 , _f - m 2 U g j 3 + ~ T T ' ^ m = l/n (3.7) Proof These results follow from ( 2 . 5 8 ) and the definition of 7(1"). We prove only (3.6), the others following in analogous fashion. We bound T ( ) below for r > ry, by 7 neglecting the term KJ , i.e. N T — Ty T{J) > Pool + ry < . Poo Inserting into ( 2 . 5 8 ) : |V¥| < I [ X H + T Y T(T - ry) dr = ^ i ^ . Therefore, V|V\P| > 0 we have XKXB + 1 . 5 7 y / g ) (XB + T y / F ) 2 X (X+ 2 - l-57y/g) (X+ T Y W ' and consequently X B < X , from the monotonicity of this function. The bound with XN comes from neglecting the yield stress terms as well as the term K . The bound N 7 using XHB comes from neglecting only the terms Pool- D Remarks: 1. The above inequalities are sharp for |V^| > 0. 2. The asymptotic behaviour as |V*| —> 00 is of most interest. For XHB, following [57], we have XHB ~ |V*| , whereas evidently XN ~ |V*| and X B ~ I V\I/|. n 50 Chapter 3. Qualitative results Figure 3.1: The function x(|V*|,if,ry,K,MOO") with parameters ry = n = 1, e = 0.55 and K = p^ = 1/2. 3. It is also relatively straightforward to provide an upper bound for x- For example, writing r(7) < 2 max{/ioo7, «7 } + TV, n implies that j(r) > 7 ("?"), where m r 7m (T) = min • T — Ty — ry l/n" 2K 2/ioo Therefore, defining XM implicitly via XMH+TY |V*| = [XM +•TY/H] r / " jL 2 Tjmir) dr, J ? y leads to x(|V*|) < X M ( | V * | ) . Further, we can see that at large r, 7m (r) T — ry 2/Joo that g X 3 2 M (XM + 1.5TY/H) 6/^*1 With the lower bounds in the above proposition, this demonstrates that x(|V*|) is linear in the limit |V*| —> oo. 51 Chapter 3. Qualitative results B e h a v i o u r of ||*|| L e m m a 1 Provided that assumptions Al, A2 & A3 hold, then the solution * of (2.63), satis- fying boundary conditions (2.67)-(2.70), lies in the space Proof Multiplying (2.63) by * and integrating over ft, using the divergence theorem: / pVtf • V t f t dn = Jn / *[pV* + S ] - i / d s t Jen / V * - S - * V - f Jsi dfi, (3.8) where u denotes the outward normal to Vt. Using (3.4), we have / *[pV* + S] t v ds = Jdii / * ( - ^ + ^ , ^ - ^ ) . i / d s Jan < V o<P \\Vp-(9<t>,9s)\\L~(dn) ) I 1*1 ds JdQ < Co||Vp-(5^fl )||ioo n)||*||ffj€ ( 8 (3.9) The first bound follows from assumptions Al & A2, using the Holder inequality. The last line follows from, f <9* L / ffidn > / |* |d0+ / ffidfi > / |*(U,i)|d£ + / =» 2 / |V*|dft > / |*|ds JO, in / |* «t|d0 o |*(0,£,i)|d£ Jdfl From proposition 2 we have that X (|V*|) > ^ | V * | . 52 (3.10) Chapter 3. Qualitative results Combining all this: , A^J|V*| dft < Co||Vp-^,5 )||icc 8n)||*||i/i + ||V-f|| ||*|| 2 € ( L2 La -- g {^}ll*llia- s {^}ll^J. 3i i F where 1*1 HI = f 1/2 7 IV*I dfl 2 JQ Since * = 0 along c/> 11*| 11,2 = 0, the seminorm is equivalent to | | * | | ||*||#i Therefore, we can find constants C\ > 0 and < ||*||/fi. d ;jJ|V*| dft di 2 < Ci||*|| i f f C > 2 i, and evidently 0, for which -C ||*||^i. H 2 (3.12) Integrating (3.12) with respect to i, we can find C3 > 0 for which C \\n%i(t) 3 < inf{f}||*||^(t) JQ t=o Z Jo We see that the integrand becomes negative if Ci l*IM*)> c 2 and consequently ||*||#i(i) is bounded for all t > 0. I(u) • and dl[u] We denote by Vj the subspace of H (Q) containing functions that satisfy boundary conditions l (2.67), (2.68), (2.72) and (2.73). The space Vj is non-empty, since for example ** G Vj, where ,T, IJ. *\ P0Ut(<P) - P , , . , ^ P~ Pin{4>) * = #(</>,i) — —- + y (4>,t)- j— TTT, PoutW ~ Pin(<P) PoutW) - Pin{<P) ,T,* Tt n i and Pin(<f>) & Pout(4>)y x out a r e / n i m (3.13) the density at the inflow and outflow, £ = 0, Z, respectively. Note that the boundary streamfunctions, & *out, must both satisfy * (l,i) = * (l,i) = Q(i), in oui 53 Chapter 3. Qualitative results *i„(0,t) = *out(0 t) = 0. ) We also denote by V/o the subspace of i / (fi) containing functions that are zero at dfl and 1 equipped with the HQ(Q), inner product. Note that L (fl) 2 but for notation we work with Vip is a Hilbert space and in fact Vj^ = Vifi. Obviously Vj is an affine space, i.e. Vi = V + V , Ifi for any £ Vf. For <J>* G V/ and u € V^o we consider thefollowingfunctional, I(u): I[u] := { hlu] + h[u] u € Vi,o, (3-14) where I [u] = f L (Vu) k k dfi, k = 1, 2, (3.15) with i l ( V u ) = L2(V«).= 21 ^ ^ d s + V^'+^-f, (3.16) -^|V** +Vu|. (3.17) Except where stated, we shall regard ** as fixed. Later, where we make comparisons between solutions, we shall explicitly denote the dependence on The subdifferential of I, di, is defined as follows: 8I[u] := {v e H \n) : I[w] > I[u] + {v, w - u) V w G V/, }, 0 where (•, •) is the inner product in Proposition 3 is The functional 0 (3.18) L (tt). 2 I(u) is strictly convex, proper and lower semi-continuous; monotone. Proof 54 dl[u] Chapter 3. Qualitative results For convexity we may follow the same steps as in [57], (proposition 2). Note that L(Vu) = Li(Vu) + L (X7u) is convex. Therefore, (see [17]: theorem 3.1, p. 66 2 and theorem 3.4, p. 74), L(Vw) is weakly lower semi-continuous. The lower semicontinuity follows from the weak lower semi-continuity. Lastly, since I[u] is convex, proper and lower semi-continuous, monotonicity of dl[u] follows from [24], (theorem 1, p. 524). 3.1.2 • Definition of 9 J I , £ 2 [ M ] and dI ,L [u] 2 2 We now consider characterization of dl[v]ij2. We consider u,w € 1 h[u] is Gateaux-differentiable. Consequently, 91^2[u] = I[ /iM 2[u], L H$(Cl). The functional and we find that > h[u] - J (w - u ) V • ( > j ^ J ( v * * + V«) + f ) dil n (3.19) for u G V] o, \/w G Vjfl: If u G Vifi, then dl ,L M 2 2 * characterized by: s (ry ' \~H (V** + VM) |V#* + Vu| (see [1]), which is set-valued for |V** + Vu| = 0. Thus (3.20) We now combine these 2 terms. P r o p o s i t i o n 4 Let I[u] = Ii[u] + I [u], 2 valued. with h[u] Gateaux-differentiable and dli[u] single Then if v\ G dKi[u] v = v\ + v G d(h[u] 2 + I [u\) & v 2 2 G dl [u). 2 Proof The sub script L means we are working with the L norm, and the first index implies e.g. I\. J 2 2 55 Chapter 3. Qualitative =>: Because results and h[u] is Gateaux-differentiable, we have v\ G dli[u) > (tu-u,-ui), Vtu e / 1 M - / 1 H i.e. because <9ii[u] is single valued. We know that v\ h[w] - h[u] + h[w) - h\u] >{w-u,v + + v G d(h[u] + /2M), 2 2 thus (3.22) GV . v) 1 (3.21) F/,0, Ifl Adding (3.21) and (3.22) we get: h[w] - I [u] >{w- u,v ) Vwe V7,. 2 Therefore, v G 2 dl [u). 2 <=: By assumption v 2 dl\{u] 2 Thus, it follows that for fixed I[w] i.e. v G dl i[u] is single valued and v\ G dl\[u}. (3.22) and (3.23) are satisfied. Hence G dl [u], (3.23) 0 2 Vw G u G V>, , 0 Therefore, because v + v G <9(/i[u] + I [u\). x 2 2 • V : Ifi > I[u] - f{w- u)V • (S[u] + f) dO, (3.24) Jn 4 « E —V • (S[u] + f). h We may now consider the steady problem for \&, which may be written as: V • (S + f) == 0, (3.25) with boundary conditions (2.67), (2.68), (2.72) and (2.73). 3.1.3 S t e a d y state 'Theorem 1 ij) = \j>* _|_ U g ) problem [Steady state problem] a n d U s i foe minimiser s There exists a unique solution * s G Vj to (3.25), where of: inf I[v}. (3.26) vev t]0 Proof This is essentially the same result as in [57] except that, due to adoption of (2.51), we are now in the Hilbert space setting. • 56 Chapter 3. Qualitative 3.1.4 results Transient problem To prove existence and uniqueness of the transient model we need to define the following operator. Let p satisfy the assumptions A1-A3, and consider the elliptic problem E[z] = -v, where E\z] = V • [pVz], in ft, z = 0on<9ft. Because for v G < p(c) pi m n < m i7 (ft), this problem has a unique solution z _1 8.7 pg. 171), and we define (3.27) is an strictly elliptic operator with bounded coefficients, thus E p ax, • E~ : l if (ft) -1 H-> HQ(£1) € H^(ft), (see [29], Th. 8.3 and Corollary as this solution, i.e. £ [?J] = z, _1 i.e. if p = constant, then i ? = A . -1 - 1 Note that because V-(S[u]+f) is the characterization of dlip.[u], this means that V-(S[u] + f) C L (ft), 2 which is compactly embedded in iJ (ft). _1 From (3.24) we can write: /(u;-u)V-'(S[u]+f)dn JQ = = [ (w-^EE^V-(S[u} + f) dft, JQ / (to - u)V • [pV(.E V • (S[u] + f))] dft, _1 JQ integrating by parts we have: f(w- JQ u)V • (S[u] + f) dft = - / V(io - u)pV(E V • (S[u] + f)) dft, -1 JQ because of the definition of E~ . l 57 Chapter 3. Qualitative results We may note that since 0 < p in < P < Pmax-, the associated seminorm defined by the inner m product, < x,x > (v,w)v = / pVv-Vwdfl, pQ is equivalent to || • Vu,weVio, Jn because of the zero boundary conditions. Note that this is just a weighted seminorm in HQ(Q) which we shall denote || • ||v" op Therefore the subdifferential of I[u] with V o seminorm and inner product is defined as: Pi dlrrx [it] := {v G V : I[w] > I[u]+ <v,w - u>Vw EVIQ}. With characterization: I[w] > I[u] + f pV(w - M ) V ( £ - V • (S[u] + f)) df2, 1 Jn i.e. v G dlfy [u] => v £ E- V l P r o p o s i t i o n 5 If u e • (S[u] D ( d l i 2 ) + f). then u G D(dlv )pfi Proof Suppose it G D ( d l i 2 ) ; which implies that 3w G dlj?[u]. Thus I[ty] > I[u] - =» /[to] > 7[u] - (w- Jn f(w- Jo u)v dtt \fw G V> ^EE^v dfl Vio G V/, )0 0 =*• J[io] > J[u] + / pV(w - u)V(E i;) dQ _1 Jn Vu; G Because of the definition of the operator E[z] we have, I[w] > I[u] + / pV(w Jn - u)Vz 58 dtt Vto G Vifl. V IFI Chapter 3. Qualitative Therefore z G dIv [u] Pi0 • (S[u] + f). E^V 3z G dlv [u] p0 3.1.5 results and because v G By definition, this implies that we can characterize dl^[u] u G D(dl), provided that u G D(dlv )- dl[u] ^ z by, z G 0. Because • pfi Existence of a solution We are now in a position to demonstrate existence of a solution. Consider the differential equation: u'(t) + A[u(t)] 3 0 t > 0, (3.28) u(0) = uo, where uo G H is given and if is a Hilbert space. ^4 = 5/ is the subgradient of I, which may be nonlinear and perhaps multi-valued. Our result follows from the application of the following result. Theorem 2 For each UQ G D(dl) there exists a unique u G C([0, co); H) with u' € function L°°(0,oo;H), such that, 1. u(Q) = u, 0 2. u(t) e D(dl) 3. -u'(t) for each t > 0, anal G di for a.e. t > 0. Proof For a proof see [24] pages 529-533. • Remember that we have two characterizations of di one with the norm || • ||^2 and one with the seminorm || • ||y . We showed that 0 u 0 G D(dI 2) L then u G 0 dI i\u] L D(dI )Vf):0 59 ^ 0. Then, by Proposition 5 we have that if Chapter 3. Qualitative results Therefore, applying Theorem 2 to this initial condition, we have that there exists a unique function u £ C([0, oo); Vjfi) with u' £ L°°(0, oo; V7,) 0 such that, 1. u(0) = u , 0 2. u(t) e 3. -u'(t) for each t D(dIv ) Pt0 £ dI , Vp 0 for a.e. 0, and t > 0. This implies that -u'(t) £ E~ V • (S[u] + l -u'(t) > f), thus £ E - ? • (S[u] + f) -E[u'(t)] 1 € V • (S[u] + f). Therefore (2.63) has a unique solution u £ C([0, oo); V/.o) and, (see earlier comments), Vr,o is 3.2 Continuity 3.2.1 Preliminary results We commence with some qualitative results for the derivative of \ with respect to some rheological properties. Behaviour of Proposition 6 J^, For a fixed gap-averaged speed, |V*|, and given £ {TY,I,TY,2\, then (3.29). Proof Now A = x(|V*|) + ry/if is given implicitly by: 60 Chapter 3. Qualitative results Because |V*| is fixed the derivative with respect to any parameter q is: -2 A 1 0 = . AH / A'H j(AH) / X J f _ , 2 f AH >- + ,i. . r '(r) dr, 7 1 /-A.H M# A J r 1 2 ^(r)dr + fAH fAH 4'' _ f / 1 ( r ) d T + _ y T y ( T ) < J r , (3.30) where the ' denotes differentiation with respect to q. Thus, — f r '(r) dr , ^ / ' . AH A' = A[ 7 V (3.31) Jry Differentiating the constitutive law with respect to ry r(j) = Pool + Kj + Ty, n at fixed r, we find: <9 1 dry Poo + 7 ' ,n—1 KUJ note that this is true for every r £ (ry,AH]. Differentiating the constitutive law with respect to r we have: fry _ 1 dr poo + Kni ~ ' n l Let / ( ) = Moo + Knj ~ , n 7 1 then A' < — ry x' + ^ < - (3.32) +4 Ty H (3-33) H (3.32) holds because r( ) > ry. Therefore 7 X'<— ry 61 '•• • (3-34) Chapter 3. Qualitative results Proposition 7 For a fixed gap-averaged speed, | V * | , and given K* e [KI,K2]> then tl K K OK Proof Differentiating the constitutive law with respect to K, r(j) = Pool + Ty, K-y + n at fixed r, we find:. ?j 1! = 8K <n, POO + Knj n note that this is true for every r e (ry,AH]. 1 Recall that, 1 dj Or poo + rcn-y™ ' -1 As in Proposition 6, let f(j) = Moo + Kn'y n \ and note that ->r. T K Thus, A' ^x' < — (3.36) K < -K +HK ^f 62 •• (3-37) Chapter 3. Qualitative Proposition 8 results For a fixed gap-averaged speed, | V * | , and given p^ € [p oo,iiMoo,2]j then < + dpoo Moo (3.38) d x Hp, Proof Differentiating the constitutive law with respect to n, T(i) = Pool + Kj + TY, n at fixed T , we find: _ * L '= dpoo 2 Moo + ^"7™ < o, 1 note that this is true for every r 6 (ry, AH]. The proof follows in an analogous way as in Proposition (7), A' ^x' < — < J L €H , 2 forC •• ( - °) 3 4 then (\A\ + \B\f- (\A\ - (A, p + WMoo L Moo Proposition 9 Let (3.39) Moo p 2 A-B) + \B\ ~ (B, P 2 B - A)) > C\A - B\ 2 (3.41) < 1 and 1 < p < 2. Proof From now on we will denote a = |A|, 6 = \B\, c = \A — B\. Because we have an inner product we can write: (A, A - B) = \a\\c\ cos a, '• 63 Chapter 3. Qualitative results where a is the angle between A and A — B. In a similar way we have: (B,B - A) = \b\\c\ cos 8, (A,B) = \a\\b\ C 0 S 7 , where 8 is the angle between B and A — B and 7 the angle between A and B. Figure 3.2: a) 7 > TT/2, b) 8 < TT/4, 7 < TT/2, C) 8 > TT/4, 7 < TT/2. Case a) 7 > TT/2, Figure 3.2a. Clearly (A, A - B) > 0, (B, B - A) > 0, then \A - B\ = (A, A - B) + (B, B - 2 A), and because (A, A — B) > 0 and (B, B — A) > 0 we have: (A,A-B) = \A\ - (A,A-B)\A\ P 2 <\A\ - (A,A-B){\A\ 2 P P 2 + \B\) - 2 P (3.42) (B,B-A) = \B\P- (B,B-A)\B\ -P<\B\ - (B,B-A)(\A\ 2 2 P 2 + \B\) -P 2 (3.43) adding (3.43) and (3.43) we have the result with C = 1. 64 Chapter 3. Qualitative results Case b) B < TT/4, 7 < TT/2, Figure 3.2b. For this case, clearly (see Figure 3.3) Figure 3.3: 3 < TT/4, 7 < TT/2 cosac > a + —, b>a+d=a+ where a = rr — a. Thus, (|i4| + |B|) -P(|i4|P- (i4, A — B) + \B\P~ (B, 2 2 2 B - A)) \A-B\ = (a + 6) - (a - cos a + 2 p p > (a + 6) - (6 " cos /? - a 2 p p cos /?) 1 1 p_1 cos /3) = (a + &)" (& - a)(p - l ) ^ cos /3 for b e [a, 6] 2 p -2 = (a - + (2 - p)6 &)(& - a)(p - l ) ^ cos/3 2 p 1-p -2 •> (2-p)6 - (6-a)(p-l)6 - cos/3 (fr-a)(2-p)(p-l) 2 > JVc, 2 p p where TY = 1/4(2 - p)(p - 1) < 1. Case c) /3 > TT/4, 7 < TT/2, Figure 3.2c. 65 2 Chapter 3. Qualitative results For this case we express cos (3 as: cos j3 COS(6J — 7) = = cos a COS 7 + sin a sin 7 ~ cos a — ^7 cos 7 + sin 0:7 ~ cos 2 a + 2 sin 57 because 7 is small. Then, (|i4| + \B\) 'P(\A\P- (A, 2 A - B ) + \B\P~ {B, 2 2 B - A)) \A-B\* (a + 6) - (& - cos/3 - a?"1 cos a|) c (a + b ) - P / V ( b P - + ( y - a?- )! cosa|; 2 p p 1 1 2 > 1 1 1 7 > - c sin 7 c Nb sin a b N > ~ 2 where ./V = 1/2. • Proposition 10 Let fields c\ and c respectively, 2 *2 <£ ^(O), constant for (2.63)-(2.65) for concentration then j ( ^ j ^ f (V*i) - ^ | p ( /or a 6e too solutions V * 2 ) ) • (*i " * ) V 2 C which does not depend on $ior Proof 66 \& . 2 d n * C||V(*i - * ) | | £ 2 2(n) , (3.44) Chapter 3. Qualitative results Now, assume that exists f i C fi such that + V * i • V ( * i - * ) > 0 and V * • V ( * - * i ) > 0 in 2 2 fl 2 + then, using Proposition 2 we have Therefore, x ( i y * i | ) _ v + x(|v* i) ^ 2 |VWi| _ v ^ ^ c | v ( ^ _ ^ 2 ) | 2 _ ( 3 4 5 ) |VW | 2 Without loss of generality assume that V * i • V ( * i —* ) < 0 in some region £l~ c fi, 2 and just for notation let us write x i = x ( | V * i | ) and x = x ( | V * | ) thus, 2 X l r V * i • V ( * i - * ) + 7=^7-7 Vtf • V ( * - * i ) = . |v*i| 1 v 12 , |v* 2 • v(*i - * ^ ^ . ( i v ^ r v f r i + |V* |- V* -V(* -*i))fJ 2|V* | 1 X2 |v*i|- v*i • v(*i - * 2|V* l-e e 2 2 2 | + 2 ) + i v * X 2 2 2 1 1_£: 2 £ 2 I V ^ I - V * ! • V ( * ! - * )) 2 for 0 < e < 2 2 r v * •v(* 2 X l 2|V*i|!- ) + | v * r v * £ e 2 - *i)) 2 2 (Ij^pj -\^^) 2 • V ( * - *i)) 2 (3-46) 1. Now we shall prove that . j v ^ ( 3 - 4 7 ) is an increasing monotone function. Let us take the derivative of (3.47) with respect to IV*|, thus what we will have to prove is x'-(l-e)T^T7 >0. ' V* v 67 Chapter 3. Qualitative results Using the fact that \ l s linear as |*| — > oo (see remarks after Proposition 2) and that, X~^|V*|, we have that in the limit |*| — > co, , X" 6/Joo because x £ C°° with respect to |V*|. Therefore, X' - (1 - e ) ^ ~ ^ ( 1 - (1 - e)) > 0 (3.48) Now, using Proposition 2 and (3.7) we can find that in the limit |*| — >0 |V*|~C7 m + 1 X , thus, x X ~ C|W| c 1 / ( m + 1 ) L—iv^i /^" " ) , 1 ~m + V 1 1 -1 then ' - (1 - e)-£- ~ C I V * ! / ^ ) - (—!—• - (1 - e) 1 X 1 (3.49) 1 which is greater than zero if £ >1- 1 m+ 1" For values of a < |V*| < 6, we have, Kmin — X ~ Kmax Bruin — X ~ B m a x , because x £ C°° with respect to |V*| and we are in a closed interval, thus x ' - ( l - e ) p ^ y >B -(l-e)K , min 68 max • (3.50) Chapter 3. Qualitative results which is greater than zero if e>l-—. max Iv Therefore, X |V#|i-e is an increasing monotone function for some range of e sufficiently close to 1. Using the assumption that V * V ( * i - * ) < 0 in 0~, we have that | V * | > | V * i | R 2 2 in Or, then we have V * i1 • V v( *1 i - * ) + T ^ - T V * , X L |V*i| X " |V* | R 2 ( X l 2 VlV^il " 1 + •V ( * - *i) > 2 2 6 + (|V*i|- V*i • V ( * ! - tt ) E |V* | 2 2 1 iv* r v*2-v(* -*i)) £ 2 2 (3.51) Using (3.5) in Proposition 2 and Proposition 9 we have, 1/ Xl , X2 1 + | V * | - 2 £ V * 1 2 •V ( * 2 6 2 - *!)) > Therefore, the result (3.44) holds. 3.2.2 v*ir v*i-v(*i-* ) £ 2 V i v ^ i i - ^ • IV*!!! - c{||^±||||p|V(*i - * )| (3.52) 2 2 • Continuity with respect to rheological and physical properties Consider the flow of a fluid with a fixed set of physical and rheological properties along the annulus, i.e. these may depend on (</>,£). Also, we are interested in the limit e —> 0, therefore the concentration is fixed in time. Suppose that one of the fluid properties varies continuously, for example from a density field p\ to a density field p , while the other properties remain 2 unchanged. We would like that our solution * i also varies continuously to *2, * i being the solution for p\ and * the solution for p . We present below a number of such continuity 2 2 results of the solution with respect to the rhelogical parameters p, ry, K and p^ and the initial conditions. For the rest of this section we will assume that the boundary function \I>* is Lipchitz continuous with respect to the rheological parameters. 69 Chapter 3. 3.2.3 Qualitative results Steady State Density continuity Note that for this case the only density dependence comes from the boundary function (3.13), The function x depends on all the rheological parameters, which we will fix for this case, but density. The variational formulation of the steady state problem is: J ^ ^ | ^ ^ ^ p ( V * * + Vu)) •V(u-w)+^(|V**+V7j|-|V**+Vu|)+f-V(u-u) dfl > 0. (3.53) Let v = wi, meaning that ui is the solution of our problem for p(a) and we denote the solution for p(c ) u = ui, then: 2 J (*jw f n ( W 2 ) 2 ) ' V ( U 1 " U 2 ) + V * 2+ V M l 1 _ 1V * 2 , ) + f 2 ' V ( M 1 ~ U 2 ) d Q ~ °" ( 3 ' M ) Now, let v = U2 and u = u\, then: J ( [W!| X l ) ( W l ) n ) ' " ~" V ( 2 W l ) + + 2 l " I V * i | ) +fi • V ( « - «i) dft > 0. (3.55) Vw 2 Note that, /^(|v*5 + v«i|)dn < /"^(|v*5-v*n + |v*i|)dn (3.56), [ T F ( | V * I + Vu2|)dft Jfi -n < / ^ ( | V ^ - V * t l + |V*2|)dft. Jn n (3.57) Substituting (3.56) and (3.56) into (3.54) and (3.55) respectively and adding the resulting variatonial inequalities, we get: + 2 / ^(|V*5 - WJ|) + (f - fi) • V(ui - ua) dfl > 0. 2 70 (3.58) Chapter 3., Qualitative results Now, le.t us keep every functional in terms of V*fe and V * j ^ , for k = 1,2. Thus, we have / (i^<^>-w ^») --» ^ ( V( 2)d /(W^">-i^fH-''-") a 2 v( , d0 fl + ( |v* l X l)(V 2 ^ ~ lv7 l' 2) X ) ( W l ) 1 ) • V ( * ~ ** 2 } d Q - ( 3 ' 5 9 ) Similarly, / (f - fj) • V ( m - ua) d f i = / (f - fi) • V ( * i - ¥ ) dfl + [ (f - fO • 2 2 2 Jn 2 Jn - *I) dft Jn (3.60) Let us consider the following functional: /.|V**+Vi;| 1 s / ds 1 2 x ( i/2) s This functional is convex, lower semi-continuous, Gateaux differentiable and ' v |V** + Vu| ' v Using Proposition 5.4 in [22] we have: F(«i)-F(u ) >' (F'(u ), *i - * ) F(u )-F(ui) > (F'(«i),* -*i) 2 2 2 2 2 where (F'( ), * - * ! ) = jf ( ^ y f f v * ) • ( * 2 - *i) d f i . Ul V 2 Therefore, /(w*'-wS-*-* ^' v v n 71 ! d , 0 ( 3 6 1 ) Chapter 3. Qualitative results T h u s from (3.58) we have / (f - f i ) • V ( * i - * ) d f i + / (f - f i ) • V ( t t $ - tfj) dfl + f 2 ^n | V ( * ^ 2 Jn + 2 /X1^ x(|V*i|) ( Jn\ W 2 W i Jn 2 dfl Jn "l^ )^~* _ x ( | V ^ 2 l ) | . y ^ . ^ ) o>Q. ) ( W l ) V ( t W 2 ) d f t (3.62) d | V * i U s i n g C a u c h y - S c h w a r z inequality, A l & A 4 , the L i p c h i t z c o n t i n u i t y of * * w i t h respect t o p, the facts that •vflV7\T/h and # i , # € - f f ^ ) . 1 2 (3-62) becomes - P2\\LI > J v Q * i - * | v f f V * 2 ' ' V ( U l " U 2 ) d " - °- ( 3 ' 6 3 ) Therefore, using the result i n P r o p o s i t i o n 10, we have (3.64) ||PI-P2||L* > C | | V ( * i - * ) | l i . 2 a Yield stress continuity A s s t a t e d i n section §2.4.2 we k n o w that x = x(l v"*|; H, Ty, K, n, Moo), thus we w i l l denote Xfc(|Vtf|) = x ( | V * | ; f f , 7 Y f c , K , n , / i o o ) = for* 1,2. W i t h o u t loss of generality, assume Ty\ > Ty , d o i n g T a y l o r e x p a n s i o n for x i ( | V * | ) a r o u n d 7y 2 2 a n d neglecting 0((TYI — TYI) ), we have: 2 x(|V*i|,ryi) =x(|V*i|,ry ) +^ x(|V*i|,r )(7yi-ry ). 2 y ; 2 (3.65) L e t t i n g v = u\ a n d u = u i n (3.53) we have: 2 / (* , jX*, V ^ a /) • V ( m - ua) +w ^ ( | Jn V | v w | 2 ( 2|) VM^ + 2 72 7m | - | V * | ) + f • V ( m - ua) dfl > 2 0. (3.66) Chapter 3. Qualitative results U s i n g t h e t r i a n g l e inequality, the right h a n d side of (3.66) is less t h a n or equal t o jf ( ^ ^ | ^ ^ p V ^ ^ •V( 2 W l -u )-f-^(|V*a-V*tH-|V*i|-|V* |)+f-V(ui-u ) 2 2 2 d f l . (3.67) N o w let v = U2 a n d u = u\ i n (3.53)and use the t r i a n g l e i n e q u a l i t y as i n (3.67), thus (^^p *l)•V(u V 2 -Ul) + ^(|V*5-V*I|-r-|V*2|-|V*l|) + f-V(U2-«l) dfl > 0. (3.68) A d d i n g (3.67) a n d (3.68) we get: + / Jn + (^2-TVl) ( | W i | _ | W 2 | [ ^ ^ ^ ( I V ^ i - V ^ I ) Jn )d n ti H dfl>0. (3.69) U s i n g C a u c h y - S c h w a r z i n e q u a l i t y a n d t h e fact t h a t ry £ L°°(fl), w i t h respect t o r y a n d the fact t h a t \T/ € H ^), the L i p c h i t z c o n t i n u i t y of * * we c a n find u p p e r b o u n d s for t h e second a n d 1 the t h i r d integrals o n the left-hand side of (3.69). t (ZV2_ZiVi)(|V¥ | _ | g , | ) ci H 1 Jn V 2 d < /"iTya-TyiUV^i-V^ldn Jn < ||TYI - T y | | 2 | | * i - < C||TVI - 7 v | | i ~ 2 2 | | f f i ' 2 f ( T 2 2 ± Z n l | v * i - V * 5 | dn < C f Jn Jn H * L |7y -Tyi|(|Ty +Ty |) 2 2 1 dVl < ||TVI - T y | | 2 | | T y i + < C||7V1 - 7V2||£a < C||7Vl ~ 7 V 2 | | i ~ . 2 i ry || 2 2 L S u b s t i t u t i n g (3.65) i n t o t h e first i n t e g r a l of the left-hand side of (3.69) we have: 73 (3-70) (3-71) Chapter 3. Qualitative results l i m - m l k 2 | | * i - *2|liji + C | | 7 v i - m | | i 2 > [ ( M i n , • - Jn V | v w i | |VW2| - * , ) df! > 0. (3.72) / The last inequality on (3.72) holds because now x ( | V * | ) has the same physical parameters and we can use (3.61). We can bound from above by: ^ | Vi* 9 W 2 | ) ( ( v* )-4 £^(v^ |V*i| 2 |V*I-*5|dfi+ f Jn I m - m l dry | V ( i t i - u ) | dft. 2 Using the result in Proposition 6, we get / m,i Jn < - ml dry | V ( U 1 - U 2 ) | d f t < C / |7V,i ;i Jn x(|v*i|) + - |V(wi - u ) | dft 2 C [ | 7 Y , i - 7 y , | x ( | V * i | ) | V ( u i - U 2 ) | dft + C / I - T V . I - 7 y | | V ( u i - « 2 ) | dft 2 )2 Jn < - ml Jn C f | m - 7 y , | x ( | V * i | ) | V ( * i - * ) | dft + C f |ry,i - T V , | X ( | V * I | ) | V ( * I - Jn + C /" K , ! - r y , Jn 2 2 2 Jn 2 | | V ( * i - * ) | dft + C / | 7 y , i - 7 y , | | V ( * I - ^ ) | d n 2 Jn 74 2 dft (3.73) Chapter 3. Qualitative results Using Holder's inequality and the fact that Ty^, i m n < 7Y,fc < 7Y,fc,max for k = recall that 1,2, the rheological properties approach a constant value as £ —• 0 and £ —* Z, we have / iTV.l-TV.alxdV^lDIV^-^ldfi < ||7V,l-7Y, || =o|| (|V*l|)|| 2||V(*l-*2)||La < C||7y,i-7Y, ||L» 2 L X L 7fi because *fc e / \TY,I Jn i7 (ft) 1 (3.74) 2 and x is a C°° function of |V*|. Similarly, 7y, |x(|V*i|)|V(*i - dft < 2 ||TV,I - ry,2|||oo||x(|V*i|)|| 2 L < C\\TY+ - 7y, ||£oo. 2 (3.75) because the Lipchitz continuity of **.. Now, see that / » ( l ' « . l ) _ 2 a ^ l i v , , |V*J - *5| dft |V* | ' |V*i| , Jn ^ X 2 ( | V * | ) , , , X2(lV*i|) , , ,, ^ V* | |V*i| ( m v ( 1 2 M 2 /V7 { T /V7 T 2 ) 2 < Assume that I ITVI C||7Vi-Ty2||Li»<C||7yi-7y2||L~. — 7V2||L°° < i C l then |TVI - ry2||L°° > Using this fact and combining ( 3 . 7 0 ) , l i m - TYallicc (3.76) > jf ( X 2 ^ 2 l ) | | m - m i l l (3.71), (3.74), (3.75) 0 and ( 3 . 7 6 ) we have: ( V * ) - ^ ^ p ( V * ! ) ) • V(*i - * ) dft > 2 2 0, using Proposition 1 0 we have the continuity result, U m - Trellioo > C||V(*i - * )||1 . 2 75 2 (3.77) Chapter 3. Qualitative results Consistency and High-shear viscosity continuity The results on continuity for consistency and high-shear viscosity can be derived basically in the same way as the one for yield stress continuity, thus we omit the proof and state the following results, - HAW - «2||ioo MOO,2||L~ > • (3-78) C7||V(*i - * ) | | i , 2 a (3.79) > c||v(*i - * )||£ , 2 a 3.2.4 Continuity for Evolution equation To prove continuity with respect to the physical parameters for the transient model will use the same techniques and results from §3.2.3. Density continuity The variational inequality for the evolution model is: + I ^-(|V#* + Vu| - |V#* + VM|) + f -V(v-u) dft>0. Jn H (3.80) Following the same technique as in §3.2.3, we let v = u\, u = u in (3.80), then letting v = u 2 2 and u = u\ and adding the resulting variational inequalities we have: J p V*2t • V(ui - u ) dfi - J piV*it • V(wi - u ) dfl 2 + 2 2 l(w<^»-i# *'») «'--» + (v V( / 2^(|V*| - W?|) + (f2 - fi) • V(«i - ua) Jn H df! > 0. (3.81) We already have a bound for the third integral in (3.81), the first and second integrals in (3.81) can be rewritten as: 76 Chapter 3. Qualitative results [ p V*2t • V(ui - u ) dfl2 Jn = 2 - V*u) • V(ui - / P2(V*2t Jn =I (P2 Jn + f pi W Jn u) 2 - Pi)V* • V(*i - * ) u 2 2 dft + / (p - pi)V*u • V(ui 2 Jn dfl- Jn I p (V* 2 u • V*2t) u 2 Jn lt u )dft 2 - tf ) 2 dfl - *I) dft f(P2- P i ) W • V(*^ - *J) dft - / p (V* - V* ) • JQ - • V(«i - u ) dfl u 2t / (P2 - Pi)V*« • V(* - * ) dft - I f p ^|V(*! - * )| dft Jn ot 2 x 2 2 z J O 2 f (P2 - Pi)Vtfit • V(tf$ - *I) dft - / p (V*n - V* ) • V(V - *I) dft (3.82) + 2 Jn JS2 2t 2 Combining (3.62) and (3.82) we have: ^(p - Pi)V* • V(*i - * ) dft - i j f p ^|V(*i - * )| dft 2 2 + + lt 2 2 2 / (f - fi) • V(*i - tt ) dft + / 2^(.|V*5 - V*I|) dft 2 2 Jn Jn w y (p - p i ) w • v(*^ - * D dft - y p (v*it - w 2 u 2 2 i ) • v(*^ - * j ) dft > => / (P2 - Pi)V*« • V(*i - * ) dft + / (f - fi) • V(*i - * ) dft 2 Jn + - 2 2 Jn / 2^(|V*^ - WJ|) dft + / (ps - pi)V* • V(*5 - *I) dft Jn lt Jn yp2(V*i -V* 0-V(^-*t)dft>X^||V(*i-*2)||| 2 t 2 (3.83) Because of condition, A.l p(c) e jL°°(ft). Having this in mind, the fact that ifr' (t) e L°°(0, k ^(fl), oo, if ^)) forfc= 1,2, the Lipschitz continuity of ** with respect to the 1 77 5 < = Chapter 3. Qualitative results physical parameters and using Holder's inequality, we have: / (p2 - Pl)V* U • V(*l - *2) dfl Jn < < ' ||pl-p2||i«»||V(*l-*2)||£2||V*lt||L2 C\\pi - P2\\L°° [(P2 - Pi)Wit • v(*5 - *i) dn < lbi-P2||i<»||v*i || 2 T < L C\\pi -P2III0C • / (P2 - Pl)V(*l, - *2,t) • V(*5 - *l) dfi < - P2|lioo||V(*l, - t t Jn [(f 2 - fi) • v(*i - * ) dn 2 < C\\P! -P2WI00 < cup! - < C||pi p ||i-liv(*i - * )lli2 2 2 Jn — p ||£oo 2 • / 2 ^ ( | v * 5 - v * i | ) d n < c||pi -P2IU00 Jn -« Combining all this we have, 11 Pi " P2IU- > * ^ | | V ( * i - * ) l ! | 2 (3.84) 2 Using GronwalPs inequality in (3.84), we have the continuity result: ||pi-P2||L-T + | | V ( * - * ) | | 2 ( 0 ) > ^ | | V ( * - * 2 ) | l i 2 1 2 i (3.85)- 1 Yield stress continuity Using the variational inequality as above, for two different yield stress fields, ry \, Ty,2, we have t / p(W - Vtfit) • V(tf 1 - tf ) dil + [ p(V* - W ) • V(*J - *3) dCl Jn Jn 2 2t + .11*1 - ry || 2 > 2 L 2t u Cjf (^^W) - ^ p ( V * i ) ) •V( ? 1 - * ) d<7 > 0, 2 (3.86) using the result in (3.77). Using the fact that ** is Lipschitz continuous with respect to TV and that *'(£) e H , we can find an upper bound for the second integral in (3.86), thus 1 78 Chapter 3. Qualitative results ( p(V*2t - V*u) • V ( * i - * ) dfi + / p(V*2t - V*it) • V(*J 2 - "il ^* ~ * ' 11 1 2)i 2l2+cllrly ~ dO T2ylli2 - (3 87) Combining the results in (3.77) and (3.87), we have the following differential inequality: liny - T I\ 2 > C^||V(tfi 2Y L (3.88) * )|||2 2 Using Gronwall's inequality in (3.88), we have the continuity result: liny - r y|| T + ||V(*i - * )||x,»(0) > C||V(*i - * )\\h 2 L2 2 2 (3-89) Consistency and High-shear viscosity continuity As in §3.2.3 the results on continuity for consistency and high-shear viscosity can be derived basically in the same way as the one for yield stress continuity, thus we omit the proof and we just state the results, - K \\ 2T+ 2 L ||V(*i - # )|M0) > C||V(*i 2 * )||£2 2 ||MOO,I-^OO,2||L2T + ||V(*I-*2)||L2(0)> C | | V ( * i - ^ H ^ . 3.3 (3.90) (3.91) Qualitative behaviour of solutions Here we derive, purely formally, a number of qualitative results for the steady and unsteady problems. These results do have a physical interpretation, which we give, but also serve as a useful test in verifying that our numerical procedure is correct. First we discuss the functions and *out, and how these are computed for a given flow. 79 Chapter 3. 3.3.1 Qualitative results E n d conditions and test concentrations The fluid concentration c manifests only in the physical properties p(c), ry(c),re(c),n(c) & Moo(c), via specified closure laws. We restrict our consideration to concentration fields c that become constant in some subdomains of fi, close to the ends £ = 0, Z, i.e. c —> 1 as £ —* 0 and c —• 0 as £ —> Z. This is physically sensible, when we study displacements that happen in the interior of fl. It follows that the physical properties p(c), ry(c),re(c)& n(c), are each constant in some neighbourhood of the ends. We denote these constant properties as follows: TY,in = 7V(1), Pin = Pout=P(0), lY,out = TV(0), Moo.tn = M o o ( l ) «in = K ( 1 ) , "in = " ( 1 ) , K t = K(0), n o „ t = n ( 0 ) , Moo.out = M o o ( 0 ) . o u The end condition ^i (<f),t) is the stream function that corresponds instantaneously to the n steady parallel flow of a fluid with properties (pi ,TY,i , n Kim " i n j Moo.in), n through an infinitely long annulus, at flow rate Q(t). To compute such a steady parallel flow, we observe that * and all flow variables are independent of £. Therefore, we have that ^ ( S , + /,) = <>, and since is a constant, this implies that 5^ is independent of 0. Now, if \I> is independent of £, we have that: <9I \> S^ = (x + Ty/tf)sgn{—}, so that the axial speed is of one sign and that A = x + TY/H is independent of <f>. We may recall that A was the absolute value of the modified pressure gradient. We now substitute from §2.4.2 the closure laws |V*| = |V*|(^4; H, ry,re,n, / i o o ) , and therefore write: f <9\I>f _^d0=/ Jo Cl(p Jo 1 Q = % (l)= n 1 |V*|(i4;ff,7y i i n ,K i f l ) n n,/ioo,in) i d^. (3.92)' The above is a nonlinear equation for A , which has a unique solution since we have seen that |V*| increases with A , when the yield stress is exceeded. We denote by A = Ai the solution n of (3.92), and now define \I>j as follows: n *in(</>) = r<t> Q\$ • / Jo -^rr t>= ci<p (i( Jo r<t> |V*|(A ;H.Ty.in, in 80 Ki ,n ,/xoo.in) n in d0. (3-93) Chapter 3. Qualitative results Similarly, * ut(0, t) is the stream function that corresponds instantaneously to the steady paralO lel flow through an infinitely long annulus of a fluid with properties (Pout,TY,out, out, out, Moo,out), K n at flow rate Q ( t ) . We then define * i in exactly the same way as we have * j , but using the OU rheological properties (ry t, K )0U o u t n , nout, Poo,out) • 3.3.2 Steady state problem In the classical setting, we characterise the steady problem as the solution * G Vj of s (3.94). 0 = V-(S + f). As discussed previously, this may be formulated more precisely as * = ** + u , where u G VQ s s s is the solution of the following minimisation: inf I[v], (3.95) v€V 0 where I[u] is defined in (3.14). Note that if we choose ** = then evidently u = 0. Equally, s from our consideration of the subgradient of I[u], we have that Vw G Vb: I[w] - I[u ] > - J(w a - u )V • [S[u,] + f) s dft. Multiplying (3.94) by \& , integrating overft,and using the divergence theorem: s 0 = f Jan * S[uJ-«/ds+ s /* V-f-[x(|V*s|)+7y/fr]|V*s| Jn S dft, (3.96) from which we see that I(u ) B = - f Jn L ( | V * | ) dft + / VsS{u } Jan 4 S s • (3.97) ds v where L (|V* |) 4 S 1 i r i v v s r v( fs s ^)l r\™s\* 1 / 2 X(|V*S|)|V* | - - J S X O ^TTj^ds > 0. (3.98) The function in (3.98) increases monotonically with |V*s| and is strictly positive for |V*s| > 0. Proposition 11 The function in (3.98) increases positive for |V*g| > 0. 81 monotonically with |V*s| and is strictly Chapter 3. Qualitative results Proof Taking the derivative of I/4(|V* |) with respect to |V* |, we have: S S L = x'(|v* |)|v* | >0, 4 s a with a minimum when |V* | = 0. S Therefore L4(|V* |) is an increasing function and is strictly positive for |V* | > 0. S S • Comparison results We consider a displacement between fluids 1 and 2. The fluid properties are specified smooth functions of the concentration: p(c), Ty(c), K(C), n(c) & Moo(c), for c G [0,1]. These functions interpolate between fluid 1 properties, at c = 1, and fluid 2 properties at c = 0. For example suppose that for tile same two base fluids we have two different fluid concentration fields, say ci and c . These may for example correspond to different stages during a displacement. The 2 two corresponding steady solutions are denoted ^s,i and *s,2We assume that the flow rate Q is identical for both flows under consideration and that the end boundary conditions on * are computed in accordance with the procedure described in §3.3.1. Since the end conditions are dependent only on the (constant) geometry and on the rheological properties of the two base fluids, it follows that the same conditions, are attained by and ty t, ou and ^s,2 at each end of the annulus. Therefore, in homogenizing the boundary conditions, we may take the same ** for each concentrationfield,ci and c . The 2 consequence is that if Vs,i = * * + 115,1 and if ^s,2 = ** + u ,2, then u i is in the test space S a> for «s , and vice versa. i2 Lemma 2 Under assumptions AI-A4, I[u ,i,ci\ 8 tions: 1. ry(ci) < ry(c ) 2 a.e. (<£,£) G fl. 82 < I[u ,2,C2\ for any of the following S situa- Chapter 3. Qualitative results 2. K(CI) < K(C ) a.e. (fag G ft. 2 3- Moo(ci) < Moo(c2) a.e. (fag) G ft. In each case above, it is assumed that all other physical properties remain equal between the two fluids. Proof The proof rests simply on the monotonicity of x(|V*|) +Ty/H of the above. For example, let us assume that if ry (c\) < Ty(c ) [X(|V*|) + TY/H]( ) [x(|V*|) + < CI with respect to any then 2 TY/H](C ). 2 Then we have: I[u ,ci] = sA inf I[v,5i] < inf I[v,c ] veVo 2 veVo = I[u ,c \. sfi 2 This result then follows, provided that x(|V*|)+ry/i? is monotone in the directions indicated, for fixed IV*! Now is given implicitly by: A = x(|V*|) + ry/H f 1 • ' *' V = AH ~A* J ^ ( R ) At fixed |V*| we differentiate with respect to any parameter q: -2A' f . . . . AH ^ / A'H j(AH) 2 ° = " 7 5 - ^ 7 ( r ) d r+ 1 1^-L r - (r)dr+^/ + <~ - J ^ AH r '(r)dr, 7 r (r) dr, 2 7 7 where the ' denotes differentiation with respect to q. The first integrand is positive, and therefore AH r sgn(A') = -sgn ^ r '(r) dr^ . 7 Differentiating the constitutive law, at fixed r, we find: 1 dj dry Moo + « n di dn di dpoo 7 n _ 1 n Moo + nnj ~ n 7 Moo + « ; n 83 <o, 7 l n_1 7 <o, < 0. Chapter 3. Qualitative results Consequently, A increases with any of these parameters. • Remark: It follows that, for increasing yield stress, consistency or high-shear viscosity, the functional - / L (|V* |)dO+ / * S[u ] •»/ds Jn Jon 4 5 s 8 increases. When the ends of the annulus consist of regions of constant physical properties, the end contributions to the above boundary integral are negligible, since = 0, and this integral simplifies to the drop in the modified pressure, along 0=1, say Ap, where Ap= [ L Jo S^idt = [ -^ L Jo °^ + pcos[3dt = p(0)-p(L)+ f pcosPdZ L Jo Therefore we have that: - /L (|V*s|,ci)dfi + gAp(c )<- / L ( | W | , c ) dfi + QAp(c ), Jn Jn 4 1 4 s 2 2 under the conditions above. The above functional may be interpreted physically as the total energy dissipation rate in the annulus. Therefore Lemma 2 states that the dissipation rate increases with any of the yield stress, consistency or high shear viscosity. This is of course physically intuitive. 3.3.3 Transient problem The main purpose of this section is to show that the solution of the transient model (2.63) goes to the solution of the pseudo-steady model (2.1) as t — > oo. First we treat the problem analytically, showing that in fact — * ||x,2(n) decays as t — > oo. Afterwards, we compute s numerically the solution for both systems and present the behaviour of ||* — \&S||L (O,)2 For simplicity we assume that the displacement in the annulus is steady and that the concentration field does not undergo any temporal changes, i.e. we have a fixed concentration field. 84 Chapter 3. Qualitative results Analytical Results For the rest of this section we assume that \&* = for the transient model. Thus, first we multiply (2.1) by (* — \T/) and integrate by parts to get S note that because = * then (* - * ) s s e H^(fl). Now we multiply (2.63) by (\& — \&) and integrate by parts, thus a + Jn[ ^ ( | V * | - |V*|) + f • V ( * - *) dft > 0. f l s (3.100) Using the fact that x(|Vtt* + V"l) (|V** + V«|) IV** + Vu| is the Gateaux derivative of We can use Proposition 5.4 in [22] and have the following result: F(*)-F(* ) s > <F'(* ),*-tf > s s where dfi. 85 Chapter 3. Qualitative results Thus, ( F ( * ) - f " ( * ) , * « - * > <0, , B because of Proposition 10. Therefore, adding (3.99) & (3.100) we have / pV^t • V(* - tf) dfi - C||V(* - * )||| > 0. s Jn s 2 (3.101) Because ^f- = 0, actually (3.101) states that -c||v(* - 9 )\\h > ~T|/ P s - *)l n 2 dn. (3.102) ' (3.103) Therefore, using GronwalPs inequality we have ||V(*-* )|| <||V(*-* )|| (0)es £a s i2 Kt where Pmin Numerical Results Here and for all the numerical results in the rest of the' thesis we consider the case when p ^ = 0. This is simpler to work with than p ^ > 0. As stated in §2.4, when Moo = 0 we are no longer in a Hilbert space setting, instead the steady model has a solution which lives in the Banach space W 1 , 1 + n(Q). The solution of the transient model (2.63), is in if (fi) for finite time. 1 On the other hand, recall that numerically we work in finite dimensional subspaces of iJ (Jl) 1 -in the case of the FEM we consider basis functions in H^^fl). In chapter 4 we use a Chebyshev spectral method, as an outline (in chapter 4 we fully explain the method), we express our 86 Chapter 3. Qualitative results solution as the sum of the first n Chebyshev polynomials which are C°° functions. Thus, both methods use finite dimensional subspaces of iif(f2), which is compactly embedded in 1 W n(Q.). 1,1+ We solve the system (2.63) using an augmented lagrangian algorithm, [26, 31], which fully describes the unyielded regions of the flow, if any. For a full description of the algorithm we refer to chapter 6. For space discretisation we use a finite element method with triangular elements and linear basis functions. For time discretisation we use the implicit Backward Euler method which comes out naturally from the augmented lagrangian formulation. As an example of the decay property of the transient solution, we compute the evolution of \I/(<^,£,£) from an initial condition, #(0, f, 0) = sin(7T(/>) sin(£7r). Figure 3.4 a) shows the decay of ||* — * || 2(Q) s L as time evolves. The rheological and physical parameters for the problem are: TV,I = 1, K\ = 2, p\ = 1, mi = 2, ryp = 1, «2 = 1, Pi = 1, m = 1, e = 0.4 and [3 = 0. Figure 3.4 b) shows the fixed concentration field, remember that we 2 are interested in the changes of the velocity field on the fast timescale. Rheological parameters are interpolated as in (2.53). Figure 3.5 shows the contours of the stream function \f/ for times t = 1, 5, 10 and the contours of the steady state solution with the theoretical results stated in (3.103). 87 Clearly this simulation agrees Chapter 3. Qualitative results Chapter 3. Qualitative results Chapter 4 Stability o f m u l t i layer flows As explained in chapter 1, for certain parameter combinations the annular displacement is ineffective. The displacement front advances faster on the wide side than on the narrow side and consequently develops a long finger. At long times, this flow will resemble a parallel multilayer flow, along the length of the annulus. It is of practical interest to understand if such flows are stable or unstable and that is the aim of this chapter. We use classical methods of hydrodynamic theory to analyze the linear stability of a parallel'multi-layer base flow. (i) In §4.1 we present the steady parallel multi-layer solutions that we consider as the base flows for our stability study. We.identify 4 different types of base flow: — Bothfluidsare yielded at the interface. — Bothfluidsare yielded at the interface but a small static channel forms in the narrow part of the annulus. — Fluid in the narrow part of the annulus is completely static. — Fluid in the wide part of the annulus is unyielded at the interface. These base solutions are governed by 9 dimensionless parameters. For the case where the narrow side is completely unyielded, we find an analytical expression to determine the minimum value of the interface position such thatfluid2 will be static and unyielded whenever the interface position is greater than this value. (ii) In §4.2 we give the general linear stability problem, which is derived by considering a perturbation about the base state. The base state and perturbation are solutions of the 90 Chapter 4. Stability of multi layer Hows interface tracking model, described in §2.5. The general stability problem will depend on 14 dimensionless parameters. (iii) In §4.3 we present stability results for the parallel flow of a singlefluidalong the annulus. We show that theflowis always linearly stable. (iv) In §4.4 we use the method of normal modes and derive the eigenvalue problem that governs stability of the general case of a parallelflowof 2 fluids. We consider the 4 types of base flow identified in §4.1. We show that linear instabilities can arise only if bothfluidsare yielded at the interface. (v) In full generality, the eigenvalue problem derived in §4.4, is governed by 14 dimensionless parameters. This makes any analytical understanding difficult. Therefore, in §4.5 we consider the simplified case of a concentric annulus,firstlyfor 2 Newtonianfluidsand secondly for 2 power lawfluids.This greatly reduces the parametric dependency and we are able to derive analytical expressions for the eigenvalues. This has use both in giving insight and as a test problem for our numerical code. •(vi) In §4.6 we analyze the full problem numerically. First we describe the method used, afterwards we present stability results for yield stress fluids. As expected, if the annulus is concentric, the results are comparable to the ones for Newtonian and power lawfluidsin an concentric annulus. In general it appears that introducing eccentricity has destabilizing effects. (vii) From our numerical results in the previous sections, we observed that all significant changes occur for long wavelengths. Thus in §4.7 we study the long wavelength asymptotics of the stability problem. For long wavelengths it is possible to derive a semi-analytical expressionforthe eigenvalues. We also present an industrial application and make some comparisons with the results found in [59]. Wefindthat interfacial instabilities are predicted to occur only in the parameter space where the model in [59] predicts an unsteady displacement, i.e. a longfingeron the wide side. 91 Chapter 4. Stability of multi layer Sows (viii) Finally, in §4.8 we investigate the effects of a spatial developing instability for a single fluid and two fluids along the annulus. It appears that any non trivial spatial perturbation admits spatially growing instabilities for the single fluid problem. For the two fluid problem we give a condition for instability in terms of two eigenvalue problems. We have not analyzed this problem further and will be considered as future work. In section §2.4.2, and for the qualitative results in the previous chapter, we have considered the limit e — > 0 in which the convective timescale is much longer than the viscous timescale. Here we wish to study the effects of transient coupling between the interface and velocity field. Therefore, we shall not assume e C l . 4.1 Parallel multi-layer flows We study the stability of steady multi-layer flows, that arise as one fluid fingers past another during an annular displacement at unit flow rate. These basicflowsare parallel to the £-axis, hence * = *(</>). We assume that the annulus is occupied by twofluids:fluid1 occupying and fluid 2 occupying <j> e {faA\- the rheological parameters offluidk as r y, rrik <j> S [0,fa) k parallel solutions of (2.83). Hence S = k The interface is denoted & K for k k 1,2 . Our solutions are steady = 0, M We denote = (5^,0), and (2.83) implies that: (Sk,4>,Sk^) ^5 = <j> = fa. (4.1) for the basicflow.By definition we have: _ dp p k cos (3 , *--Tz--sr-> . - Sh (4 2) and with (2.89), this implies that the axial pressure gradient is continuous at the interface: dp |2 ~d~z = 0. Therefore the fluid layers are acted on by a constant pressure gradient, modified by an axial static pressure that is different in each section due to the density jump. Since S ^ is independent k of <j> and £, (by definition, we look for solutions independent of £), following [59] we may write: S W = A, S, 2 (p 92 = A-b, (4.3) Chapter 4. Stability of multi layer Hows where b is the buoyancy parameter, given by: P2~ Pi (4.4) =— cos 8. St* Note that typically, in a cementing scenario wherefingeringoccurs, 6 is negative: the heavier , b H cement channels past the more viscous drilling mud that is left behind on the narrow side. The constant A represents the modified pressure gradient in the axial direction, withinfluid1, and must be found as part of the base solution. In order tofindA we use (2.85), with Q = 1 (for a fixed unit flow rate), i = *( 1 ) = / r<t>i d$> Jo r 1 d<t> (4.5) dfa d4> fc=2. d<p + i 4 fc= where, assuming thefluidsare Herschel-Bulkley fluids, 0, = sgn(A) { fc=i \A\ < H mi+2 \A\ ~ T /H) T /H, 1IY MI+1 IX /c (mi+2) 11 1 \A\ + \A\ > mi + ny/H, \A-b\< 5* - H {\A-b\ m2+2 sgn(A - b) { fc=2 (4.6) TI,Y/H 7i£)"i2+i T ,Y/H, 2 (4.7) T2, /HY K^{m 2 + 2)\A-b\ i \A - b\ > T ,Y/H. 2 see [59]. Equation (4.5) is a nonlinear equation for A . It is straightforward to show thatthe integrals on the right-hand side increase strictly monotonically with A and that (4.5) has a unique solution A = A(fa), depending only on the interface position, on b, and on the rheological parameters of the twofluids.Given each A(fa), we can construct * from (4.6) & (4.7), by integrating with respect to cj> from <f> = 0. We shall denote the basic solution, constructed as above, by ^k,o, (or simply *o if only a single fluid is considered). In general tykflifi) is continuous in <> /, with a jump in thefirstderivative (hence the axial velocity) at fa. We present some characteristic basic flows below. 93 Chapter 4. Stability of multi layer flows From (4.6) k (4.7) we can see that the basicflowdepends on 9 parameters, ry k, «fc, t mis, f° r k = 1, 2, b, e and </>j. The buoyancy itself depends on 4 parameters, p , St* and /?. Although k these do not appear in the baseflowdescription, these parameters will appear individually in the stability problem that we derive later. In addition, this stability problem will contain a wavenumber and the time scale ratio e. Thus, in total our stability problem has a 14 parameter dependencies and this is the main difficulty of extracting useful information from the results. For simplicity we will consider the buoyancy as one parameter, and for the rest of the thesis we fix the Stokes number as St* = 0.1. Figure 4.1 shows the base parallelflowof two generalised Newtonianfluidswith rheological parameters, Ty,\ = 2, «i = 2, m\ = 2, ry,2 = 1, « = 1, wi2 = 1, 6 = —9.5, the eccentricity is 2 0.4 and the position of the interface, denoted as fa, is at <f> = 0.5. The inclination of the well is j3 = 7T/10. AS we can see, bothfluidsare fully yielded, the pressure gradient and buoyancy forces are large enough to over come both Ty\ and Ty - The stream function is continuous and 2 a discontinuity on the axial velocity is observed at the interface. Figure 4.2 shows the basic parallelflowof a couple of fluids with rheological parameters, ry,i = 0.5, Ki = 1, mi = 2, Ty,2 = 0.6, K = 1, m,2 = 1, b = —9.5, the eccentricity is 0.5 and fa = 0.5. 2 The inclination of the well is (3 = TT/10. The pressure gradient and buoyancy forces are not sufficiently large to fully over come Ty along the narrow part of the annulus and a small mud 2 channel forms. Figure 4.3 shows a worst case scenario. Here the rheological parameters are, TY,I = 1, «i = 1, mi = 2, T Y j 2 = 10, K2 = 1, m 2 = 1, b = —28.5, the eccentricity is 0.8 and fa = 0.5. The inclination of the well is (3 = TT/10. Again, we increase the eccentricity of the annulus to 0.8. and clearly the pressure gradient and buoyancy forces are not large enough to overcome TY2 and fluid 2 if completely static. Figure 4.4 is of purely mathematical interest, to our knowledge no experimental work has reported a case like this. The rheological parameters are,Ty:i = 7, K\ = 7, mi = 2, Ty,2 = 0.2, K.2 = 1, m = 1, 6 = 0,the eccentricity is 0.2 and fa = 0.5. For this case we have a fully vertical 2 94 Chapter 4. Stability of multi layer flows annulus. Note that if we increase the eccentricity, the region where fluid 1 is flowing gets wider and the plug zone will break. Figure 4.1: Multilayer base flow, both fluids yielded at the interface. Rhelogical and physical parameters are: Ty,i = 2, « i = 2, p\ = 2, m\ = 2, T y = 1, «2 =T, pi = 1, m = 1, e = 0.4, 2 i 2 4>i = 0.5 and /3 = TT'/IO. I 1 Figure 4.2: Multilayer base flow, both fluids yielded at the interface, mud channel in the narrowside of the annulus. Rhelogical and physical parameters are: ry i = 0.5, K\ = 1, p\ = 2, m i = 2, 7 Y = 0.6, «2 = 1, P2 = 1, m = 1, e = 0.5, ^ = 0.5 and [3 = TT/10. t I2 2 95 Chapter 4. Stability of multi layer flows Figure 4.3: Multilayer base flow, fluid 2 unyielded at the interface. Rhelogical and physical parameters are: iy,i = 1, «I = 1, p \ = 4, mi = 2, ry,2 — 10, K = 1, P2 — 1, m = 1, e = 0.8, 2 2 & = 0.5 and B = TT/10. Figure 4.4: Multilayer base flow, fluid 1 unyielded at the interface. Rhelogical and physical parameters are: Ty i = 7, K\ = 7, p\ = 1, m\ = 2, T y = 0.2, K = 1, p = 1, m = 1, e = 0.2, <?> !i = 0.5 and /3 = o'. t ) 2 96 2 2 2 Chapter 4. Stability of multi layer Hows 4.1.1 Static mud channels For certain combinations of dimensionless parameters, it is possible for one of the fluids to be stationary, as we can see in Fig. 4.3. Physically, this occurs when the modified pressure gradient A , (or A — b in fluid 2), is not large enough to overcome the yield stress of the fluid. Since the modified pressure gradient is constant in each fluid layer and since H(<f)) decreases, yield stress fluids become static at the largest values of <§> in each layer, i.e. where the gap is smallest. It is therefore possible to have flows for which the fluid 1 layer or the fluid 2 layer is either partly or fully static. From an industrial perspective the flows within this category that have most interest are those in whichfluid2 is static, either partly or fully, i.e. there is a static mud channel. Let us consider a basic multi-layerflowof twofluidswith specified physical properties, computed as described above, and how the type of flow may vary for different choices of fa. Since the minimal annular gap is at <f> = 1, where H = 1 — e, there can only be a static "mud" channel influid2 if there is a static layer as fa — + 1. The condition for this to happen is that: \A(fa = l ) - b \ < ^ - e . (4.8) Suppose that (4.8) is satisfied. We now ask what is the maximal azimuthal width of static mud channel that can exist. If the mud is static, then <9* f<Pi 1 = /o io <t>\ d (4.9) dfa fc=i i.e. all the unit flow rate mustflowthroughfluid1, <f> G (0,fa\. Evidently as fa — > 0, if the mud remains static, then A(fa) —> oo since the unitflowrate is forced through an increasingly narrow azimuthal gap. It follows that for somefa,there will be a solution to ^ - h \ = wty - (4 10) which defines the smallest interface position for which the mud layer remains fully static. We denote this interface position by fa = fa, i , m n and the maximal azimuthal width of mud channel is therefore 1 - fa in- Let us now explore the parametric dependency of (4.8) and (4.10). m 97 Chapter 4. Stability of multi layer flows First, suppose (4.8) holds, so that there is a static layer. For fa > fa,min, the pressure gradient A(fa) is determined implicitly by (4.9). Evidently, sgn(A) = 1 and for AH > r^y we can rewrite (4.6) as ( A H _ \ d* n,Y\ m i + (AH 1 \ i,Y \n, MI T Y „, k=l mi + 2 which at each <f> depends only on e, mi, Ajr\y and + mi 1 +1 (4.11) ,'AH\ J Jl,Y B\ = T\y/«i- Therefore, we see that (4.9) implies A(fa) /(e,Bi,mi). T\,Y Consequently, dividing (4.10) by r^y, we have T2,Y f(e,Bi,mi) H(fa^ in) m (4.12) n,Y i.e. depends on (e,B\, mi, At, Ay), where A = h ^ T"2,y Av (4.13) n,Y T\,Y are the buoyancy ratio and yield stress ratio, respectively. The condition that (4.8) holds is now: (4.14) Ay > (1 - e) /(e,Bi,mi) n,y i.e. we have static mud channels if and only if the yield stress ratio exceeds a critical value that is dependent on (e,Bi,mi, A(,). Figures 4.5-4.9 show the effects of eccentricity on the magnitude of <&, in- As expected, if we m increase the eccentricity of the annulus a mud channel will be left during the displacement in the narrow side of it. Rheological parameters play an important role as well, clearly if ry2 > ry,i, ; i.e. Ay > 1, the pressure gradient necessary to overcome the yield stress of bothfluidsat the interface will not suffice to overcome the yield stress offluid2 in the narrow side of the annulus. Having as a result a decrease of the domain wherefa^m= 1- The same phenomena will occur if either n > «i and p > Pi2 2 98 Chapter 4. Stability of multi layer Hows Figure 4.5: Contour plot of <^,mm> eccentricity vs. yield stress of fluid one, buoyancy parameter equal zero. Rhelogical and physical parameters are: KI = 0.5, p\ = 1, m\ = 1, Typ = 0.5, K = 0.4, p2 = 1, rn = 1 and 3 = 0. The zero level curve corresponds,to the equality of, condition (4.14), i.e. <^ j = 1, the rest of the level curves correspond to different values of &,min when condition (4.14) is fulfilled. 2 2 )m n 6 <1 • ' A -(l-e)|l{e,B m )-Aj>0 Static layer A -(1-e)|l(e,B,m)-A l=0 No static tayer v ] 0.2 ] v B 0.4 I1 I 0.6 Figure 4.6: Contour plot of 4>i, in, eccentricity vs. yield stress of fluid one, buoyancy parameter negative. Rhelogical and physical parameters are: K\ = 0.5, p\ = 1.1, m\ = 1, ry^ = 0.5, K2 = 0.4, p2 = 1, m,2 = 1 and 8 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. i ^ m i n = 1, the rest of the level curves correspond to different values of r^.min when condition (4.14) is fulfilled. m 99 Chapter 4. Stability of multi layer Hows 0.9 0.8 0.7 0.6 0 0.5 0.4 0.3 0.2 0.1 0.2 0.4 0.6 0.8 1 \ Figure 4.7: Contour plot of <fo m, eccentricity vs. yield stress of fluid one, buoyancy parameter positive. Rhelogical and physical parameters are: KI = 0.5, p\ = 1, m\ = 1, ry^ = 0.5, K = 0.4, P2 = 1.1, m = 1 and 8 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. <?!>i i = 1, the rest of the level curves correspond to different values of <pi min when condition (4.14) is fulfilled. im 2 2 im n } Figure 4.8: Contour plot of 4>i i , eccentricity vs. yield stress of fluid one, equal consistency. Rhelogical and physical parameters are: «i = 0.5, p\ = 1, m\ = 1, TY,2 = 0.5, « 2 = 0.5, p = 1, m 2 — 1 and 8 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. Si = 1 the rest of the level curves correspond to different values of ^^min when condition (4'.14) is fulfilled. <m n 2 m i n 100 Chapter 4. Stability of multi layer flows 0.9 0.8 0.7 0.6 " 0.5 0.4 0.3 0.2 0.1 0.2 0.4 0.6 Ay 0.8 1 Figure 4.9: Contour plot offa nin,eccentricity vs. yield stress of fluid one, consistency of fluid 1 less than consistency of fluid 2. Rhelogical and physical parameters are: K\ = 0.3, p\ = 1, mi = 1, Ty-,2 = 0.5, K2 = 0.4, p2 = 1, m = 1 and [3 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e.fa,mm= 1, the rest of the level curves correspond to different values offa minwhen condition (4.14) is fulfilled. lT 2 } 4.1.2 Summary of results Solutions of a basic flow consists of solving a 1-D nonlinear equation, followd by integration of an algebraic relation tofind\&. This computations depend on 9 different parameters, rygt, Kfe, rrik, for k = 1, 2, b, e and fa A case of industrial interest is when a static layer forms (see §4.1.1). Such flows occur only for interfaces fa > fa, m, and for certain combinations of m parameters. We have shown that fa m depends only on 6 parameters, ry^, for k = 1, 2, «i, t mi, b and e. Therefore, provided that m fa > fa,nu n the flows depend only on a reduced subset of the dimensionless parameters. 4.2 Stability of parallel multi-layer flows For the rest of this chapter we will work with the interface tracking model, defined in §2.5. We consider the stability of any of the parallel multi-layerflowsin §4.1. We denote the transient 101 Chapter 4. Stability of multi layer flows interface by, = 0. Our system of equations can be written as: F = <p — fa(£,t) V • [S + p V*fc,] = 0, fc fc t ~d!T ~dqV~dJ ' = *i(0,^,t) = 0, *2(U,t) = 1, %lfc= = 2 (Sfc + PfeV*fc,) • n|JL 2 t 0, (4.15) mfl , k ^4> = fa(gt) • (4.16) (4.17) (4-18), at <£ = <&(£,£), (4.19) = -6(l,tan/3sin7r<fo) • n, at <£ = <&(£,*), (4.20) We consider linear perturbation of these flows, assume normal mode expansions and transform to an eigenvalue problem, in the usual way. We denote the base flow solution from §4.1, by *fc — ^kfl(4>) & Sfc = (Sfc^o,0), with interface position fa = fap, (i and recall that Sk,<j>,o constant in each domain. For < 5 < 1 we assume the expansions: *fc = ^k,o + 5^k,i+S ^>k,2 + - , 2 Sfc = (£fc,,/,,o,0) + JSfcj + <5Sfc2 + ... 2 ) fa = fa,o + Sh + ..., substitute into (4.15-4.20) and expand in powers of 5. V* fc = (*fc,o,* + ^ M , * + -.^fc,U + -) (- ) 4 21 |V*| = [(*,o,* + <J*M,* + - ) + (***,U + -) ]' 2 2 fc ~ I*fc,o,*l+^gn(*fc,o,^*fc,i,^+0(^). (4.22) When the fluid is yielded we find: S ,i = [x (|*fc,o,*l)^,i,*. (Xfc(l**,o,*|) + W ) T | ^ T ] ! fc 1 M fc (- ) 4 ) 23 but otherwise Sfc is indeterminate: We see that the stability problem depends on 14 dimensionless parameters. First of all the base flows depend on the following 9 dimensionless parameters: T y ^ , K , mjt, forfc= 1, 2, b, e and k 102 is Chapter 4. Stability of multi layer Hows <f>i. Now the buoyancy parameter b is defined in terms of the two densities, p , the inclination k [3 and the Stokes number,, St. The densities and the inclination appear individually in.the above stability problem and therefore the full stability problem depends on 12 dimensionless parameters. Additionally, rather than solving the initial value problem above, we will take a normal mode expansion of the linear perturbation. This will introduce a 13th dimensionless parameter, the wavenumber in the axial direction. Finally, in the kinematic equation for the interface motion appears a 14th parameter e. Since this general situation is difficult to analyse, we start with some of the simplified scenarios. First we consider the stability of a single fluid along the annulus. 4.3 Single fluid problem We commence with theflowof a singlefluidalong the annulus, and assume that the base flow of the singlefluidis everywhere yielded, although the arguments following will also apply if the baseflowis unyielded in a static layer on the narrow side of the annulus. The parameter dependence reduces greatly for this problem, now we only have 5 parameters to take into account: ry, K, m, e and the wavenumber a. The linear stability equation for the single fluid is: p(A*!) = V • x'(l*o,*|)*i,*, X(l*0,*|) + T7 t (4.24) with boundary conditions: * (0) = 0 = *i(l). x We suppose a normal mode solution for * i : tf i ~ /(</>)e *- . i(c st) Substituting \I>i into (4.24) we have: isp(D 2 - a )} = D[Df '(\* M 2 X 0 / 103 X(f*0,*|)+ jy ~ a2 (4.25) Chapter 4. Stability of multi layer Hows where D = Multiplying (4.25) by /*, the complex conjugate of / , and integrating (by parts) with respect to <j>, we obtain: is[J +a h] = 2 x (4.26) J + aI, 2 2 2 where /%l/l d* 2 -fi = (4-27) k _ /V* *¥V* * ( l d JO Ji = h = (4.28, l*O,0| / p|£/| d</, (4.29) I (4-30) 2 Jo Jo \Df\ 2 X '^ all of which are real positive definite. If we suppose that the normal mode is temporal, i.e. a > 0, then clearly s = is/ is imaginary and the growth rate s/ < 0, so that the flow is linearly stable. Therefore, temporal instability can only occur for the multi-fluid problem, i.e. due to the presence of an interface and the jump conditions found there. In Figure 4.10 we show the spectrum of problem (4.25) with rheological parameters, Ty,\ = 1, K\ = 1, mi = 1, pi = 0.5, on the left-hand-side and p\ = 1, on the right-hand-side. We show our results in the following order, first we plot the spectrum as a —> 0, from (4.25) we may simply derive the analytical bound for Sj, the upper and lower analytical bounds are defined by: — maxx' < Si < —minx', (4-31) The argument of x' is |*o,<^|(^) d the max and niin are taken over <f> £ [0,1]. It has been a n shown that x' > 0, thus it is clear that the imaginary part of the temporal normal mode is always negative for a small . Second, the pictures in the middle of Figure 4.10 show the maximum of the imaginary part of the eigenvalues for varying wave number a. 104 Chapter 4. Stability of multi layer flows Figure 4.10: Spectrum of the temporal normal mode with lower and upper bounds as in (4.31) & (4.32), figures on the leftrhand-side show a fluid with rheological parameters: left-hand-side, 7 y = 1, K\ = 1; p\ = 0.5, mi = 1, e = 0.4 and 8 = 0 on the right-hand-side we only change p\ to 1. (1 105 Chapter 4. Stability of multi layer Sows Finally, we show the spectrum as a — > oo. Here from (4.25) we can derive: max(x + TY/H) min|^| - + Ty/H) max|*' | ' min( < X 1 - [ 0 ' ' as we know, X + TY/H has been defined as the absolute value of the modified pressure gradient, thus a positive quantity. 4.4 Two-fluid eigenvalue problems We have seen above that any single fluid parallel flow is linearly stable to spatially periodic temporal perturbations. We ask therefore, whether the presence of an interface between two parallel streams can linearly destabilise this flow. The linearised equation for V • [Sfc,i + AbVtyfc,i„i] = 0, is inftfc, (4.33) with Sfc,! defined by (4.23) when yielded. The kinematic equation and boundary conditions are: —h + #fc,u + ^kfldh = 0, at <p = <&o , (4.34) *i,i(<U,t) = 0, (4.35) * ,i(U,i) (4.36) = 0, 2 The jump conditions (4.19) k. (4.20) are linearised, first about the basic flow and secondly onto the basic flow interface position. Thus, at 0(6), condition (4.19) becomes (*fc,i + /i*fc,o,*)|fcI 2 = 0, at 0 = ^,o, (4.37) and we note that it is the derivative of this quantity with respect to £ that appears in (4.34), i.e. it is irrelevant which fluid is considered for the kinematic condition. Expanding (Sfc + Pk^^k,t) about (j) = 4>ifi w e have: [Sfc + PfcV*fc,] • n ~ t Sk,o,(j> + Sh^Sk,o,4> + 5Sk,i,(i, +Pfc[*fc,o,tft + Sh^kfiMt + ^k,i,4>t\ s at <j> = (I,tan/3sin7r0j) • n ~ 1 — Sh^ tan 8sinnfafi 106 / fop n, Chapter 4. Stability of multi layer flows Noting that Sk,o,<t> 1 S independent of <f> and ^k,o independent of t, at 0(5), condition (4.20) becomes (S'fc,!,^ + Pfc*fc,i,^i)lfci = Wig tan/3sin 71-^0, at <f> = = 4>ifi, (4.38) We consider first the case where both fluids are yielded and mobile throughout the annulus and subsequently those cases when a fluid layer contains static regions. Mobile fluids In this case Sk,i is defined by (4.23) and the stability problem, consisting of (4.33-4.38), is posed on the entire fluid domain. We assume a normal mode solution: i,fc and substitute into (4.33)-(4.38) to give: T\,Y is (D -c?)h 2 Pl (4.39) £[X2(l*wl)I?/2]-a isp {D - 2 2 2 o?)f 2 (4.40) 0, (4.41) 0, (4.42) and at H sh + a(fk + */c,o,0^o) = 0, (4.43) (fto*o,M + /fc)li = °> (4.44) 0 [(Xk(\Vo,k,<j>\)-ispk)Dfk]\lJl = (4.45) iabho tan (3 sin ir<f> ifi In general, it is necessary to solve the eigenvalue problem (4.39)-(4.45) computationally. We discuss results for this in §4.6. 107 Chapter 4. Stability of multi layer Hows We note that for temporal stability we may consider a^O, since if a = 0 then (4.43) implies that either s = 0 (in which case /fc = 0 and the perturbation is trivial), or instead ho = 0. For = 0, following the procedure in §4.3, we may multiply by /£, integrate across each fluid layer, sum the resulting expressions and show that s = isj with sj < 0. The interfacial terms cancel by virtue of (4.44) & (4.45), which simplify considerably. Static mud channel problems: 4>i > 4>i,min In this case fluid 2, on the narrow side, is completely unyielded and hence static. Therefore, for (f> >fawe have 32,0 < H(4>iY By choosing 5 sufficiently small we may ensure that |S ,o + 5S ,i| < 2 2 T2.Y H(faY so that that fluid 2 remains unyielded. In this case, the interface does not deform under a linear perturbation and we have again a single fluid problem: xi(l*wl) + ^ l*i,o,*l ! D[xi(|*i,o,*IWi]-a 2 > r i fi, = is {D -a )h z (4.46) 6(0,^,0) /i(0) = 0, /i(&,o) ' = 0, l Pl ' (4.47) (4.48) As in §4.3, we can multiply through by / j \ integrate over ( 0 , ^ , 0 ) , and establish that these flows are linearly stable to temporal perturbations, but unstable to any non-trivial spatial perturbation. Partly static fluids In this case we first suppose that fluid 1 is fully yielded and fluid 2 is part yielded and part unyielded. We discuss afterwards the other possibilities. Therefore, we have a yield surface at say 4>2,Y where <t>2,Y = 02,y,O + &/>2,Y,l(£, *)> 108 Chapter 4. Stability of multi layer flown with fayfi the yield surface of fluid 2 in the base multi-layer flow and perturbation of the yield surface. Note that fa < fay < 1 and fa < is the linear 5fay \(£ ,t) t > fa, i m n As 4> —>fat fluid2 becomes unyielded across the annular gap and static, i.e. Y 1X7*21-* 0 <j>-^fat .. as Y Expanding each component to 0(5) and linearising onto * 2 , u + <A2,y,i(^i)*2,o,^ = * ,u fayfl, 0, = °' 2 we have that: at <f> = fi at (A = </>,y,o. Equation (4.50) leads to the boundary condition, / = 0 at 2 (4.50) 2 <j> = (4.49) fay , fayfl, for the normal mode. This condition is the same as (4.42). Equation (4.49) simply describes evolution of the yield surface once the stream function perturbation is known. As in other linear stability analyses of yield stressfluids,[27, 28], the yield surface perturbation plays an essentially passive role, in that the same boundary conditions are satisfied with or without the yield surface. On the other hand, its worth noting that the yield stress affects the stability problem; via the basic flow and (here) via definition of the yielded region. Our eigenvalue problem is identical to (4.39)-(4.45), but replacing (4.40) & (4.42) with : T2,y ^2 | ^ 2 0 4> \) £>[X2(l*wl)Dh\ - o? —— ^l*2,0,*lr ^ - l / 2 , 4> G = isp2(D'-a*)h (fafi, Hfayfi) fayo) = 0, (4.51) (4.52) i.e. only the domain of thefluid2 problem is changed. The other situation to consider is that in whichfluid1 is partly unyielded. Since H(fa) decreases, the criterion |Si| < T \ y / H will be met in some layer <p G (<Ai,y,o,fa)for the basicflow.It follows that |Si | < Tiy/H at the interface, for the baseflow,and hence that for sufficiently small 5 109 Chapter 4. Stability of multi layer flows the interface does not yield under the perturbation. The stability problem in fluid 1 becomes: Xi(l*i,o,*l) + • i V111M/1 n^,i) -t- U[xi(|* |)I>/i]-a 2 W R h, l*i,o,*l = ispi{D -o?)h 2 0€(O,0i,y, ) (4.53) O with boundary conditions (4.41) and /I(0I,Y,O) (4-54) = O. For fluid 2, we have either (4.40) or (4.51), depending on whether the base flow of fluid 2 is fully or partially yielded. In either case,fc{4>i)= 0, since the interface is not deformed, and either (4.42) or (4.52) is satisfied at <f> = 1 or <j> = <j>2,Yfi, respectively. Neither of these cases is of any interest. For both, we may follow §4.3, multiply the fluid k equation by /£, integrate across the domains and establish that the flows are linearly stable. Summary To summarise, the only base flows that may be temporally unstable are: (i) two fully yielded fluid layers; (ii) fluid 1 fully yielded and fluid 2 partly yielded (i.e. a static channel on the narrow side, wholly contained within fluid 2). Stability of the other baseflowsessentially reduces to the stability of one or more independent singlefluidlayers. 4.5 Analytical simplifications Analytical progress is possible only in limited cases. One of these cases is when the annulus is concentric and consequently the base velocity is constant in each fluid layer. The simplest case is when there is no yield stress. 4.5.1 Newtonian fluids in a concentric annulus For two Newtonianfluidsin a concentric annulus (H(<p) = 1), we find: .6 1 + (1-</>i,o)~— *o^ = • % * W = * 4>i,0 + (1 - <Pi,0)— K2 110 «i t b 0»,OQ— <Pi,0 + '(1 -V ' <Pi,0) — K-2 (" ) 4 55 Chapter 4. Stability of multi layer flows The two fluid layers are mobile and in fluid k, Xk and |V\I/| are related linearly by: Xk 3/«fc|V*|. Therefore, (4.39) & (4.40) become (3 -Vis)CD -a )/i 0, <f> e (3«2 - ip2s)(D 0, <i> e (fa,o, 1) 2 2 M - 2 We note that if — (ink ipks) 0 and these modes are stable. Therefore, for (3«fc — ipks) ^ 0, and interestingly observe that the eigenvalue = 0, instability we may assume that o?)h then sj < (0,fa,o) s does not appear in the field equations, but only the jump conditions. This degeneracy is peculiar to the Newtonian case. The solutions that satisfy fi{fa) =Ai for constants 2 2 ^ 0, write A\ 1 1 sinh afao 2 we find + (l-fafi) fi = 2 1 «2 3/t «2 sinha(^i,o - 1) for constants «l 6 C 2 + fa a e = hrjCit ^4 = hoC & (4.44) Q Ci (4.56) 2 Ci & C , and divide (4.43)-(4.45) by h . From (4.43) 2 are: f (fa) = A sinha(0-1), sinha0, We may assume that ho A\ & A . (4.41) & (4.42) + (1 - & , o ) fa a e 0i,O 2 fi «2 Finally, substituting into (4.45), we find the following quadratic relation for s: i 3/t tan B sin nfa o (4.57) 2 •cotha^n 1 ac Kl + 0i,O + (1 - fa,o) K 2 • Pi Z S — 3K 2 K 2 b Kl « • cotha(fafi — 1) a £ 0i,O 2 + 3K (1 - 02,0 1 — is 2 0i, ) O K 2 111 P2_ 3« 2 Chapter 4._ Stability of multi layer flows Figure 4.12: Maximum imaginary part of s vs. a . Rheological and physical parameters are: « i = 1) Pi = 1) P2 = 1 e = 0.0, <fii = 0-5 and (3 = 0. 112 Chapter 4. Stability of multi layer flows Figure 4.14: Maximum imaginary part of s vs. a. Rheological and physical parameters are: = l,e = 0.0, fa = 0.55 and f3 = 0. a) p = 1, P2 = 0.75. b) pi = 1, p = 1.25. K l 2 x 113 Chapter 4. Stability of multi layer flows Figures 4.11 and 4.12 show the maximum imaginary part of s as a function of the wave length a. As we can see, buoyancy is a key parameter for instability to arise. In the absence of any density difference, it appears that viscosity differences are not strong enough to drive an unstable behaviour. Thus, it seems that we need the velocity difference at the interface to be sufficiently large for instabilities to appear. Note the difference between the results found by Raghavan & Marsden, [62], and ours. They showed that for a horizontal multi-layer parallel flow in a porous media in the presence of density difference, if the heavier fluid was in top of the lighter one, the flow was unstable for all viscosity ratios, if the case is reversed, the lighter on top of the heavier, they claim that the flow will be unstable only for large viscosity ratios (the more viscous on top of the less viscous). For the vertical parallel flow, buoyant forces play a key role in defining the jump in the tangential velocity. Figures 4.13 and 4.14 show the effect of the interface position on the stability of the flow. Clearly, effects of this position are not of great importance, for a concentric flow. Figure 4.15 shows the effect of inclination of the annulus. As expected, in the presence of 114 Chapter 4. Stability of multi layer flows buoyancy, theflowmay be unstable if the heavierfluidis on top of the lighter one, i.e. p\ > p . 2 If we have positive buoyancy, i.e. p\ < p , then the inclination has stabilizing effects. 2 Another case for analytical progress is when there is no yield stress but the power law index n is less than one. Again we will assume that the annulus is concentric and that the base velocity is constant in eachfluidlayer. 4.5.2 Power-Law fluids in a concentric annulus Because the base velocity is constant in eachfluidlayer, we find: \A — b\ \A\ mi m2 where A is the modified pressure gradient defined in (2.57) and b is the buoyancy. We can find A{4>i) from: I A\mi 14 — h\ m2 In general (4.58) has to be solved numerically. The two fluid layers are mobile and in fluid k, from the definition of Xk, a n d having no yield stress we see that: Xi{\*o,*,i\) = \A\ X2(|tto,* |) = \A - Using (4.59) we can find an explicit expression for x'k ,,. T ., < (mi + 2) Xi(l*o,*,il) |*o,*,i| = «r(mi + 2) l^l" 11 terms of A , K2™ (m + 2) .. 1 and m X2(|*0A2l) l*0A2l 1 115 (4.59) b\. l2 2 2 < (^ + 2) 2 = 2 \A-b\^-i • y • > Chapter 4. Stability of multi layer flows Substituting (4.60) & (4.61), into (4.39) & (4.40), we get _ DVi /^*>-wry\ q2 V D h-a — , = , ,, 0 e(0 (0) , (4 2) 2 -M -w 1+) 1 mi = 0, m 0G(<Pi,o,l) (4-63) with boundary and jump conditions: /i(0) = 0 (4.64) / (1) = 0 (4.65) 2 (4.66) • L * • , / m i K \ ( m + 2) 1 /K n2 2 ( m + 2) \ 2 (4.67) and the kinematic condition + a (/ /l i++ /iosgn(A) osgn( — H + jTOl ' ' , , )=0 at 0 = &, (4.68) / (0) = C /i sinhA (</>-l), (4.69) 0 0 Solving (4.62) & (4.63) with (4.64) & (4.65) respectively we have: / (0) = Ci/i sinhAi \ 1 o r 2 2 0 2 for constants C\ & C and 2 A 2 A 2 = = ^ f^(mi+^-tl^-VtS a 2 / / " (^2+2)-i]A-6r-y 2 t 2 116 2 S Chapter 4. Stability of multi layer Sows Assuming that ho ^ 0, and using (4.68) we can find an expression for C\, C l _ = ( 1 sg*(A)\A\ ' s sinhAi^j mi \ae K MI (mi +2))' In a similar way, using the first jump condition, (4.66), and the definition of C\, we canfindan expression for C , 2 C 2 s = -r sinhA 0j 2 sgn(A-b)\A-b\ * m K {m \ae n2 2 + 2) 2 Using the definitions of C\ & C and (4.67) wefinda non-linear equation for s, 2 «7 (mi+2) Ai tanh Ai<; A tanh A<j l l iabho tan 8 sin ITfao mil^l™!K ( m + 2) m A|m -1 2 n a 2 • pis 1 2 2 P2S 2 2 s ae sgn(A)\A\* K ( + 2) MI MI s_ _ sgn(A- b)\A K% {m c7e~ 2 2 -b\ * m + 2) (4.70) 10 20 30 40 50 40 (a) (b) Figure 4.16: Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1.5, K = 1, e = 0.0 and 8 = 0, mi = 1.5, m = 2. a) fa = 0.5. b) fa = 0.6. For b < 0, p\ = 1 and p = 0.75, for b = 0, p\ = 1 and p = 1 and for b > 0, p\ = 1 and p = 1-25. 2 2 2 2 117 2 Chapter 4. Stability of multi layer Bows 0.4r 0.2- 11 , , , , , 0 10 20 30 40 50 a Figure 4.17: Maximum imaginary part of s vs. a. Rheological and physical parameters are: K\ = 1.5, K2 = 1, e = 0.0, fa = 0.5 and (3 = TT/8, mi = 1.5, m = 2. For b < 0, p\ = 1 and P2 = 0.75, for b = 0, p\ = 1 and p = 1 and for 6 > 0, pi = 1 and p = 1.25. 2 2 2 Figure 4.16 shows the maximum imaginary part of the spectrum as a function of the wave number a. We have considered two interface positions, as in the Newtonian case, fa does not play an important role in the stability of the flow. For Newtonian fluids, we observed that if fluid 1 was more viscous than fluid 2, but fluid 2 heavier than fluid 1, i.e. buoyancy effects dominating viscous effects, the flow was unstable. We have the same phenomena for power law fluids and asforthe Newtonian case, it seems that in the absence of density difference the viscous effects are not strong enough to drive an unstable behaviour. Finally, Figure 4.17 shows the destabilizing effects of inclination if heavier fluid 1 is on top of lighter fluid 2. Even if fluid 1 is more viscous and with lower power law index than fluid 2, stable cases from previous figures, instability arises. In this type of flow, buoyancy and inclination are the key factors for unstable behaviour. . Now we can proceed with the investigation of the full problem, i.e. yield stress fluids in an' inclined eccentric annulus. 118 Chapter 4. Stability of multi layer flows 4.6 Numerical Results In this section we present numerical resultsforthe eigenvalue problem, D [ x i ( l * w I W i ] - « is (D -a )h 2 2 2 Pl (4.71) £[x (l*2,o>l)Af ]-a 2 isp (D' 2 - cO/ 2 2 2 (4.72) and at 0, (4.73) 0, (4.74) Pi,Q- H sh 0 0, (4.75) (Ao*o,M + Jfc)li = 0, (4.76) + a(fk + ^k,o,4,ho) = k=2 [(Xk(\^o,k,<t>\) ~ ispk)Df ] k k = l = iabh 0 (4.77) tan j3sinirfa, 0 4.6.1 Numerical Method In this section we give an outline of the method used to solve system (4.15)-(4.20). We use a spectral method, which is a global method that uses the fully discretized stability operator that later on can be supplied to a matrix eigenvalue solver which gives us the spectrum. We have chosen to use Chebyshev expansions to discretize our problem. The use of Chebyshev o polynomials, especially in bounded domains, have been very effective and accurate for the discretization of initial and boundary value problems, see [66]. Chevyshev polynomials Chebyshev polynomials may be defined in several ways, e.g. 119 Chapter 4. Stability of multi layer flows • In terms of trigonometric functions: T {y) = cos(n n arccos(j/)). • As solutions of the following Sturm-Liouville problem, (v^Ar.(y)) dy V " v -^=r (y)-o. + B dy-"""J • As a recurrence relation, T (y) = l, 0 Ti(y) = y, ... ,T (y) = 2yT (y) n+1 n - T _i(«). n For numerical reasons, we will use the trigonometric representation which is the most practical. The method is summarized as follows: • Approximate /i(y) and J2{y) using a Chebyshev expansion, N fk(y) = J2a T {y), n n and evaluate the Chebyshev polynomials at the extrema of the iV-th Chebyshev polynomial, given as yj = cos ^ which are known as the Gauss-Lobatto points. Note that y € [—1,1] and that <j> e [0,1] so we use the following change of variables so we can work with the Gauss-Lobatto points, . ... Oil) for 0 6 [0,(pi], 2 -(-fa + <Pv)y + (4>i + 4>Y) for (p G (<pi,l], where 4>i is the position of the interface between the two fluids and <py is the position of the yield surface, if any. Observe that if fluid 2 is completely yielded, then <py = 1. After discretizing in this way, we have as a final result the following eigenvalue problem in discrete form, Aa = 120 cBa. Chapter 4. Stability of multi layer flows We chose to implement the boundary, jump and kinematic conditions on the first, last and next-to last row for fi(<fi) matrices and on the first, second and last row for f {4>) 2 matrices. In the case of boundary conditions, we implement them only on thefirstand last row of B. These same rows in matrix A are chosen to be a complex multiple of the corresponding rows in B. In this .way, the resulting spurious modes associated with the boundary conditions can be mapped to an arbitrary location in the complex plane. Code validation >l "0 1 10 1 1 20 30 a 1 40 1 50 -0.2 0 1 1 10 (a) ' 20 a ' 30 1 40 1 50 (b) Figure 4.18: Maximum imaginary part of s vs. a. Inclined annulus. Rheological and physical parameters are: «i = 1.5, K = 1, m\ = 1, m — 1, b = —1 and 3 = 0. a) Exact solution, b) Numerical solution. 2 2 We validate our code with-the exact value of the maximum imaginary part of the spectrum, for 2 Newtonian fluids, which is given by (4.57). The absolute value of the error values at each a is presented in Table 4.1 and the numerical and analytical values are compared in Figure (4.18). Following, we present our numerical results for the full problem. 121 Chapter 4. Stability of multi layer flows a maxsj exact maxsj numerical Error 0.001 -1.327407374138356E-8 -2.151373954828488E-8 8.239665806901317E-9 0.501 -3.284347465794034E-3 -3.284348381546957E-3 9.157529227109273E-10 1.001 -1.256769763479928E-2 -1.256769663542340E-2 9.993758789267337E-10 1.501 -2.638888993138258E-2 -2.638889188681512E-2 1.955432531175472E-9 2.001 -4.280735057269110E-2 -4.280735611097798E-2 5.538286877715404E-9 2.501 -5.993596650918038E-2 -5.993596650695535E-2 2.225025719226892E-12 3.001 -7.634083070570413E-2 -7.634083068837617E-2 1.732795851250302E-11 3.501 -9.117814039039307E-2 -9.117814032189228E-2 6.850078837494777E-11 4.001 -1.041107832784735E-1 -1.041107832084426E-1 7.003088386969836E-11 4.501 -1.151352471298047E-1 -1.151352470885577E-1 4.124704744423724E-11 6 -1.387140792107665E-1 -1.387140792299231E-1 1.915664848972654E-11 8 3-1.560598543973128E-1 -1.560598543900940E-1 7.218753372839615E-12 10 -1.653734816055713E-1 -1.653734815152703E-1 9.030104441976050E-11 12 -1.708225692269705E-1 -1.708225692174500E-1 9.520467747492489E-12 14 3-1.742498623405992E-1 -1.742498624152761E-1 7.467687579421067E-11 16 -1.765335018738152E-1 -1.765335016295237E-1 2.442915036926507E-10 18 -1.781267104457814E-1 -1.781267104777813E-1 3.199993048319527E-11 20 -1.792802857595362E-1 -1.792802855209202E-1 2.386159325684645E-10 22 -1.801413773516178E-1 -1.801413771880445E-1 25 -1.810736771505876E-1 -1.810736770116384E-1 1.389492687575711E-10 28 -1.817271639851303E-1 -1.817271639827821E-1 2.348121697082206E-12 31 -1.822025801255159E-1 -1.822025804615425E-1 3.360265676821683E-10 34 -1.825590619823075E-1 -1.825590618802357E-1 1.020717121935633E-10 37 -1.828331366144432E-1 -1.828331366346426E-1 2.019936995445448E-11 40 -1.830483364567075E-1 -1.830483481665764E-1 1.635732937987910E-10 1.170986885234981E-8. Table 4.1: Code validation. Absolute value of the error between analytical and numerical results. 122 Chapter 4. Stability of multi layer flows 4.6.2 Yield stress fluids Figure 4.19 shows the effects of yield stress and interface position on theflowin a concentric annulus. As expected, because bothfluidsare yielded at the interface, yield stress does not appear to have a large influence on the stability of theflow,and nor does interface position. Only viscous and buoyant forces govern the unstable behaviour. The results are very similar to those of the Newtonianfluidsin Figures (4.11)-(4.14). Figure 4.20a) shows similar behaviour as for Newtonianfluids,compare with Figure (4.15). As above, bothfluidsare yielded at the interface. Figure 4.20b) shows the effects of inclination of the annulus. As in the Newtonian case, if the heavier fluidl is on top of the lighterfluid2 the flow will be unstable. (a) (b) Figure 4.19: Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1.5, K2 = 1, e = 0.0 and 8 = 0, mi = 1, mi = 1, ry^ = 0.9 and Ty i = 0.7. t = 0.6. For b < 0, p\ = 1 and p = 0.75, for b = 0, p\ = 1 and and p = 1.25. fa 2 2 123 pi a) fa = 0.5. = 1 and for b > b) 0, p\ = 1 Chapter 4. Stability of multi layer Hows (a) (b) ; Figure 4.20: Maximum imaginary part of s vs. a. Rheological and physical parameters are: a) K\ = 1.5,re = 1, e = 0.0, fa = 0.5 and 8 = 0', mi = 1, m = 1, ryj = 0.7 and ry = 0.9. b) «i = 1.5, K2 = 1, e = 0.0, & = 0.5, /? = ir/8, mi = 1, m = 1, ry,i = 0.7 and ry '= 0.9., mi = 2, m = 1.5. For b < 0, p\ = 1 and p = 0.75, for 6 = 0, p\ = 1 and p = 1 and for b > 0, pi = 1 and p = 1.25. 2 2 ]2 2 2 2 i2 2 2 Following this line of reasoning, if we have a concentric annulus and both fluids are yielded at the interface, changes in power law index will lead us to qualitatively similar results as for power law fluids, as examined in §4.5. Next we will study the effects of eccentricity in the stability of the flow. Figure 4.21 shows the effects of eccentricity on the flow. As we can see, eccentricity has destabilizing effects if 6 < 0 and stabilizing effects if b > 0. This may be explained as follow. We are considering flows for which both fluids are yielded at the interface. Now, when we increase the eccentricity, we need a higher pressure gradient to overcome the yield stress of fluid 2, note that we assume that fluid 2 is on the narrow side of the annulus. If fluid 2 is less viscous but heavier than fluid 1, the velocity difference between both fluids will decrease due to the fact that the heavier fluid is flowing in the narrow part. Having the opposite situation, the less viscous and lighter fluid in the narrow side, the velocity difference will increase and instabilities will arise. 124 Chapter 4. Stability of multi layer flows 0.15r 0.5. (c) (d) Figure 4.21: Maximum imaginary part of s vs. a. Different positions of the interface between the fluids. Rheological and physical parameters are: K\ = 1.5, K = 1, e = 0.0 and 3 — 0, 2 mi = 1, m 2 = 1, 7y,i = 0.9 and 7y i 2 = 0.7. a) e = 0.2. b) e = 0.4. c) e = 0.6. d) e = 0.8. For 6 < 0, pi = 1 and p = 0.75, for b = 0, pi = 1 and p = 1 and for b > 0, pi = 1 and p = 1.25. 2 2 125 2 Chapter 4. Stability of multi layer flows -6 1 0 .§i -0.5 1 1 1 0.5 I i (a) r- 1| , , , 1 2.5, (b) (c) Figure 4.22: Spectrum and eigenfunctions for system (4.39)-(4.40), a = 0.5. Rheological and physical parameters are: «i = 1.5, K = 1, mi = 1, m = 1, -Ty,\ = 0.9, 7y,2 = 0.7, negative buoyancy, e = 0.4 and 8 = 0. a) Spectrum of the eigenvalue problem, b) Eigenfunctions of the interfacial mode, solid"line shows the real part of f\ and f , dashed line shows the imaginary part of fi and f and dashed thick line shows the magnitude of fi and f . c) Eigenfunctions of an interior mode, solid line shows the real part of fi and f , dashed line shows the imaginary part of fi and f and dashed thick line shows the magnitude of fi and f . 2 2 2 2 2 2 2 2 126 Chapter 4. Stability of multi layer flows Lastly in Figure 4.22a) we show the spectrum of the problem with rheological and physical parameters K\ = 1.5, K 2 = 1, mi = 1, m 2 = 1, Ty i t = 0.9, ry^ = 0.7, negative buoyancy, e = 0.4 and 8 = 0 for wave length a = 0.5. Clearly, we have only one unstable mode, the interfacial one. We proved in §4.3 that the problem of a single fluid is always stable, thus if we are away from the interface all the modes are negative. In Figures 4.22 b) and c) we show the eigenfunctions for the interfacial mode and an interior mode respectively, note the big jump at the interface for the interfacial eigenfunction this may be due to the large difference in the base velocities when instability is present. Another interesting point that it is worth to mention is that all significant changes occur for long wave lengths. This is also in agreement with our model derivation. Recall that our model neglects all the behaviour at 0(6*), thus a » 1 will fall into this regime and our model loses validity. Therefore we will concentrate in the limit a —> 0 by studying the long wave length asymptotic behaviour of our problem. 4.7 Long wavelength asymptotics First, we consider the flow of two Newtonian fluids in a vertical or near vertical concentric annulus. 4.7.1 Newtonian fluids We are interested in the behaviour of the solution of (4.57) for small a, and we have seen that s— > 0 as a —> 0. Thus we assume that, s ~ as\ + a S2 + ••• 2 Let us use the series expansion for coth(aj) around zero, substitute all these into (4.57) and collect first order terms, then 127 Chapter 4. Stability of multi layer flows 0(1) problem kfi - —(<i>i,0 - 1) I K2 (<t>ifi - 1) - ) Sl + e ( *0,1,* — \ K2 *0,2,tf<fo,0 )= 0 Thus, S l = e ^*o,i,*^(l " &,o) + *o,2,*&,o) r. v 4>i,0 H '- \ Kl (-78) 4 (1 - <Pi,o) ) «2 / which is a real value. To determine the sign of Si in general is hard, since (4.78) depends on 6 parameters, «i, n , pi, P2, St* and 4>i$- Now, we continue to the 0(a) problem. 2 O(a) problem Collecting 0(a) terms, we have the following equation for s , 2 % \ 4>iflJ e\<pifl-l • 2 ( P2 1 Pi 1 \3K2 e(cj>ifi -1) \3K2<Pi,0 . 3/t e<fo, 2 0 3«2 (<pi,0 - 1) 4- «~—tanpsin7r</>jo, 6K2 Therefore, P2 , S 2 ISi I Si I = 1 0iO ~\3K ^ , 02) ^ , 0 tan/?sin7r^,o *i,o + £(l-&,o) ' ( 79) Clearly, s is pure imaginary and in the presence of negative buoyancy, i.e. fluid 1 heavier and 2 in top of fluid 2, and some inclination, i.e. 8 ^ 0 the flow may be unstable. On the other hand, 128 Chapter 4. Stability of multi layer Hows if we consider a vertical parallel flow, then the stability curves will be determined by the sign of s\. For example, let us assume for the moment that si > 0, thus Substituting (4.78) and after some algebra, we have e<pifl (— — (Vox* \Pl «2 - *o,2,*) + (*o,2,^ - ^.l,*)) > 0, (4.80) / the inequality is reversed if s\ < 0. This will lead us to 2 cases, 1.) Case 1, (*o,2,* " *o,i,*) > 0, thus b < «2 (4.81) P2 2.) Case 2, (*0,2,* - *0,1,*) < 0, thus Kl - «2 > ^ > «2 jj (4.82) P2 If Si > 0, the cases are analogous with the inequalities reversed. From these expressions we can see that it is quite difficult to determine the sign of si in a general way, the 5 parameter dependency makes the analysis very difficult. In spite of this difficulty, we 129 Chapter 4. Stability of multi layer Bows can infer that the jump in the base velocity is the principal physical cause for the interface to be unstable, i.e. if the jump is sufficiently large, instabilities arise. This translates as follows, the jump in the base flow is governed by viscosity difference, buoyancy and interface position^ thus, both viscosity and density ratio play a decisive role in the instability of theflow.It is worth to mention that this phenomena is not comparable to Kelvin-Helmholtz instabilities, see e.g. [18]. This classical instability involves two inviscidfluidsone in top of the other, here only inertia, which we have neglected, governs the instability. In [33] the authors describe what they call Kelvin-Helmholtz-Darcy instability. They consider parallel flow in a horizontal Hele-Shaw cell, one of thefluidsis a gas, the other is a liquid. They report that the instability is only governed by inertia and that the viscous effects govern the wave velocity and the growth rate only. Our problem concerns the parallelflowof two viscousfluidsin a vertical Hele-Shaw cell and even though the nonlinear inertial terms are absent, instabilities arise and are governed by density and viscosity differences. The only possible comparison between the Kelvin-Helmholtz instability and the Kelvin-Helmholtz-Darcy instability with the one found in the previous work is that without difference in velocity theflowis neutrally stable. When we move to the problem of a near vertical Hele-Shaw cell, we see that the viscosity ratio no longer plays an important role on the stability of theflowif the heavierfluidis on top of the lighter one, b < 0 case. This behaviour is seen in (4.79). If b < 0 the last term of the right hand side of (4.79) will be always positive. Numerical example Figure 4.23 shows the contours of the imaginary part of s with respect to p and /t . We have 2 2 2 fixed K\ = 1, pi = 1, e = 0.0, fa = 0.5 and [3 = 0. This is in full agreement with our analytic results. One thing is worth to mention is that when we have the isodensity case, i.e. b = 0, the flow will be always stable. This is easily seen, if = 0 and for example, from (4.55) we have, «2 130 (*o,2,</> — *o,i,4>) > 0, Chapter 4. Stability of multi layer flows and for instability to arise, from Case 1 we have, 2 ! < £ = 1. Pi K2 Therefore the flow will always be stable. • \ \\ \ \ o c c •/ / \ \V ^ " / V? / / / / / \\'8 \ Figure 4.23: Contour plot of s . Rheological and physical parameters are: K\ = 1, p\ = 1, e = 0.0, fa = 0.5 and 8 = 0. 2 4.7.2 Y i e l d stress fluids Now, let us consider the asymptotic behaviour of the full problem when a — > 0. Assume that, /fe ~ fk,o+afk,i + o?fk,i s ~ a s i + a s + ... h ~ + - 2 + och\ + a h 2 2 + substituting this expansion into (4.39)-(4.45) and collecting the first order terms we have: 131 Chapter 4. Stability of multi layer flows 0(1) problem D(x' Dfkfl) 0 = k (4.83) 0 = 0,1 /fc,o = 0 and at 6 (4.84) = H —jsiho + (f fl + ^kfirfho) k (/fc,o + *fc,o>o)li x'kDfkfilt 0 .(4.85) = 0 (4.86) = (4.87) = 0. Solving (4.83) we have: /!,O(0) (4.88) = Jo /2,O(0) Xi (4.89) /2, (0i)+ / ^d4> J<t>t Xi = O where C\ / ,o(0i) = 2 / — r 0 + ^o*fe,o>l2d Applying the boundary condition at 0 = 1, we canfindan expression for Ci, Ci where /i(0O 1 (4.90) + / (l) 2 [<>> (0) = / - d 0 1 J (0) = / - ^ 0 . L 2 JO Xl J(pi 132 A2 (4.91) Chapter 4. Stability of multi layer Hows From (4.85) we get clearly, si G R. 0(a) problem D( ' Dfk,i) = X k at 2 kfi 0,1 = 0 fk,i and (4.93) PkisiD f (4.94) 4> = fa U — (s ho + sihi) + (fk,i+Vk,o,<l> i) h 2 = (/jt,i + **,o,^i)li = (x'k fk,i-isiPkDfk,o)\i 0 (- ) 0 (- ) = D 4 4 ibfro tan/?sin71-^0• 95 96 (4.97) Integrating (4.93) once, we have: x' Dfk,i-PkisiDf k = C,. kfi 2 k using (4.97) we get, C ,2 = C ,i + ibho tan dsinn 2 2 fao Integrating (4.98) and using (4.96) we can find the solutions to (4.93), 133 .(4.98) Chapter 4. Stability of multi layer flows / ,i(0) 2 where C\ is defined in = h M - h i M l k M +j + i s i K ^ ) d ^ ( ' °) 4 10 (4.90). Using the boundary condition at 4> = 1 we can find an expression for C , i 2 C 2,i = i j / J. \ i T t -\ \ -C\isi h(<t>i) + h(l) where I\ & I2 are defined in (4.91) Using the kinematic condition g(&, ) --—s 0 2 T ) , \ , -t6fe tan^sxn7r0i,o hWi) + h(l) Tm , Ii{4>i) +1 {±) o r v (4.101) 2 and (4.95), (4.99), (4.90) and (4.101) we have ., ... , Ji(^)/ (1) , ,„ , , =,6tan/Wfro + A * + 2 f i W + A f a ( 1 ) M l T 2 J {fa)I {l) - J {l)h{4>i) , . 2 / ^ ) + I2(l)) ' (4'103) x w an expression for s , 2 ( 2 l ( clearly, s = s ,/. 2 2 Note the first term of the right hand side of (4.103), as in previous cases, if b < fluid 1 is on top of lighter fluid 2, the buoyancy has destabilizing effects. 0, i.e. heavier Turning to the second term of the fight hand side of (4.103), we see that again, the jump in the velocity plays an important role in the stability of the flow. But things get more complicated, note that to determine the sign of s2, we have to investigate the effects of Fortunately, from §4.6.2 we know 11 parameters on (4.103). that the position of the interface does not play an important role on the stability of the flow. We also proved that eccentricity has destabilizing effects on 134 Chapter 4. Stability of multi layer Hows the flow as well as inclination. Therefore we will concentrate in investigating the effects of the 8 rheological and physical parameters on the stability of theflowwhen we are in the case of a mildly eccentric vertical annulus. Figure 4.24 shows the level curves of the maximum imaginary part of s for rheological para2 meters, pi = 1, pi = 1, m\ = 1, m = 1, Ty \ = 1, Ty 2 = 1. The neutral stability level curve, 2 } t Ki = K implies that iffci> k , i.e. the displacingfluidis less viscous than the displaced one, 2 2 theflowwill be stable. In Figure 4.25 a) and Figure 4.25 b) we have decreased and increased, respectively, the buoyancy parameter. When the displacingfluid,fluid1, is denser than the displaced fluid,fluid2, fluid 2 is pushed up the annulus by buoyancy. This increases the jump in the velocity and a wider unstable region develops. Increasing the buoyancy parameter has the opposite effect, now the displacedfluidis denser than the displacing one and the jump in the velocity decreases, resulting in a thinner unstable region. Figure 4.26 shows the effects of the difference in power law index. Figure 4.26 a), m\ = 1.2 and m2 = l, can be explained as follows, fluid 1 now is more shear-thinning thanfluid2, resulting in a shift of the level curve for n\ = K to the stable region. The opposite effect happens when 2 mi = 1 and m = 1.2, Figure 4.26 b). Now the level curve for K\ = K has been shifted to the 2 2 unstable region, clear effect offluid2 being more shear-thinning thanfluid1. Figure 4.27 shows the effects of yield stress on the stability of the flow, for Figure 4.27 a)fluid1 has a higher yield stress thanfluid2, TY,I = 1.5 and 7y i 2 = 1. Because of this, a higher pressure gradient is required to overcome 7 ^ 1 for bothfluidsto be yielded at the interface. This could be seen asfluid1 being "more viscous", thus the level curve K\ = K is shifted to the unstable 2 region. If now 7 y = 1-5 and Ty \ = 1, Figure 4.27 b), we see thatfluid2 is "more viscous" and )2 t the level curve «i-= K is shifted to the stable region. 2 To the knowledge of the author, this is a phenomena never reported in the literature. For this kind offlow,even in the absence of inertia, theflowcan be unstable. All rheological parameters 135 Chapter 4. Stability of multi layer Hows p l a y a role o n the s t a b i l i t y of the flow, s o m e t h i n g not seen i n N e w t o n i a n a n d power l a w fluids. F o r y i e l d stress fluids, even w i t h o u t the presence of buoyant forces the flow c a n be u n s t a b l e . F o l l o w i n g the results found i n [59] we c a n make some c o m p a r i s o n o f our contour plots a n d u n s t a b l e regions for the cases w h e n the displacement is unstable. 3| F i g u r e 4.24: C o n t o u r plot of s , 2 m i = 1, m 2 = 1, 7V,i « i vs. K . 2 R h e o l o g i c a l a n d p h y s i c a l parameters are: p\ = 1, = 1, 7 v , 2 = 1, P2 = 1, e = 0.2, fc = 0.5 a n d (3 = 0. 136 Chapter 4. Stability of multi layer Sows Figure 4.25: Contour plot of s , «i vs. K . Rheological and physical parameters are: mi = 1, 2 77i2 = 1, TY,I = 1, Ty,2 = 1, e = 0.2, 2 fa = 0.5 and 8 = 0. a)pi = 1.1 and p 2 = 1. b)pi = 1 and P2 = l . l 1 K (a) (b) Figure 4.26: Contour plot of s , KI VS. K . Rheological and physical parameters are: p\ = 1, 2 TYi = 1, ry i2 and m 2 = 1, p 2 = 1.1, e = 0.2, 2 fa = 0.5 and /3 = 0. = 1.2, 137 a) mi = 1.2 and m 2 = 1. b) mi = 1 Chapter 4. Stability of multi layer Hows r 1 2.5 1 ? / u 0.01' 1.5 J / S 1.5 (a) 2 2.5 (b) Figure 4.27: Contour plot of s , KI VS. n , buoyancy parameter equal zero. Rheological and physical parameters are: K\ = 1, p\ = 1, K = 1, p = 1, e = 0.2, fa = 0.5 and 8 = 0. a) 7y i = 1.5 and Ty, = 1. b) Ty,i = 1 and 7y2 = 1.5. 2 2 2 ; 2 4.7.3 2 } Industrial application In this section we make some comparisons with the results found in [59]. In [59] the authors use a lubrication model to predict, for a given set of parameters, whether a good displacement takes place or not. If they have an unsteady displacement, i.e. the interface front moves faster in the wide part than in the narrow part, a long finger develops along the annulus. This situation is when our assumption of parallel flow is likely to occur. Furthermore, they investigate when a viscous fingering instability will occur during the displacement. They showed that viscous fingering happens only in the unsteady displacement regime. Figure 4.28 (from [59]), shows when a unsteady displacement will occur, U meaning unsteady displacement but no viscous fingering, F meaning unsteady displacement and viscous fingering, S meaning neither unsteady nor viscous fingering. Rheological and physical parameters are, ry i = 1.0, KI = 1.0, mi = 1, 77i2 = 2, b = —0.5, 8 = 0. It is clear from this picture that as the authors increase e, 7V,2 and n , i.e. viscosity of fluid 2 increasing, an increase of the domain of 2 138 Chapter 4. Stability of multi layer Hows unsteadiness a n d i n s t a b i l i t y is observed. I n F i g u r e 4.29 we show the level curves of the m a x i m u m i m a g i n a r y p a r t of s 2 as a f u n c t i o n fo TY,2 a n d e. W e have the same fixed rheological parameters as i n F i g u r e 4.28 w i t h fa = 0.5 a n d 8 = 0. W e capture the same behavior as i n F i g u r e 4.28, as fluid 2 becomes m o r e viscous, the unstable d o m a i n increases. Surprisingly, our unstable d o m a i n occurs o n l y as a subset of the u n s t e a d y a n d fingering d o m a i n i n F i g u r e 4.28. T h i s means t h a t i f a l o n g finger develops d u r i n g the displacement, for c e r t a i n parameter ranges the parallel flow m a y be unstable. Therefore, apart from h a v i n g a n u n s t e a d y displacement, m i x i n g m a y occur at the interface of the finger due t o the i n s t a b i l i t y t h a t we s t u d y here. It is not clear w h y there is a p p a r e n t l y no interfacial i n s t a b i l i t y for p a r a m e t e r regimes for w h i c h the a n n u l a r displacement w o u l d be steady. Note t h a t for such parameters, a base parallel flow is a perfectly well-defined s o l u t i o n of the two fluid problem. 139 Chapter 4. Stability of multi layer flows Figure 4.28: Example unsteady displacement and local instability/viscousfingeringregimes: effects of changing ry,2 and e for different K2- (a) K, = 0.5; (b) K2 = 1; (c) K2 = 1.5. Fixed parameters are: Ty,i = 1.0, K\ = 1.0, mi = 1, m.2 = 2, 6 = -0.5, /3 = 0. F, local fingering instability ; U, unsteady, but no local instability; S, neither unsteady. From Pelipenko & Frigaard (2004). 2 140 Chapter 4. Stability of multi layer flows 0.2 1 0.3 0.4 0.5 0.6 0.7 0.8 0.9 (a) • , \ • \ \ s J.0I \ \ r •\ 5.1 s ^ U i 0.2 0.3 0.4 (b) I 1 0.5 e 0.6 l_ 0.7 0.8 0.9 (c) Figure 4.29: Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = K = 2 4.8 1, 7V,i = 1, mi = 1, 77i2 = 2, 6 = —0.5, '= 0.5 and (3 = 0. a) K 2 = 0.5. b) n 2 = 1. c) 1.5. Spatial stability problem Finally, we present a collection of results concerning linear spatial instabilities. Although incomplete, the results look interesting and potentially of practical value. 141 Chapter 4. Stability of multi layer flows 4.8.1 Single fluid spatial instability Let us consider,forthe singlefluidproblem, a spatially developing instability: s > 0 and a = aij + ia/. In this case, (4.25) constitutes an eigenvalue problemfora & C. We find: 2 2 2 R s Jxh + hh = Ci >0 s I +1 s[Jih - J /i] = c sI +I 2 -(a ) = a j-a 2 R (a )/ = lajaji 2 2 2 2 2 2 2 2 2 (4.104) (4.105) Solving for a h a/, we have: R ICx + y/C'f + ai = ± a C 2 (4.106) 2- = ± R (4.107) It appears that any non-trivial complex eigenvalue a and eigenfunction / must admit spatially 2 growing instabilities. 4.8.2 Twofluidsspatial instability In an analogous way as in the previous section, let us assume that: s > 0 and a = a + ia/. R We assume small s and study the asymptotic behavior of the problem. Thus, let /fc a ~ //c,0 + s//c,l + S / / c , 2 + - - 2 i i 2 h ~ ao + sa\ + s a ... ~ ho + sh\ + s h + 2 2 2 substituting this expansion into (4.39)-(4.45) and collecting the first order terms we have: 142 Chapter 4. Stability of multi layer flows 0(1) problem - a o A , o = 0, D(x' Dfkfi) k 0 A,o = at 0 = 0,1, ao(/fc,o + *fc,o,4>M = 0, (/fc,o + **,o,*'>o)li = 0, Xfc-D/fc.olf = iao/iotan/^sinyT^o. (4.108) System (4.108) has solutions of the form, fkfi{4>) = fkfi, {4>) + fkfi,h(4>), P where D( ' Df , ) x k kfi P =0 /fc,o,p = 0 at 0 = 0,1, /fe,o,p = -*fe,o,0^o- Therefore, (4.109) /l,0,p = -/l0*2,0, /2(0) 1 (4.110) '/ (0i)' 2 where = / Jo 1 "7^ Z" 1 1 J Xi 2(0) = / v d0. J<t> X 143 2 Chapter 4. Stability of multi layer, hows Now, fko,h is a n eigen-function of D(x'k fkfi,h) = <4fk,o,h, D with OQ being the eigen-value. Therefore, from Sturm-Liouville theory, oft is real and less than zero. Thus, , ctR =0 Our solution, fk,o{4>) = fkfi,p{4>) ho f*W7S + *i,o,*7^) and has to satisfy the second jump condition, then + /fe,o,/i(<^) + 2Df2,o,h X 0. aj ~ x'iDf , lfi h = -aih t nf3sm7r4>i, (4.111) 0 & equation (4.111) defines ho. Therefore, condition for instability is that A = — a - is an eigen-value of both: 2 D( [Df)=\f x =o /(O) = D(X2Df) fifo) = A/ = /(l) = 0. For 2 Newtonian fluids at a fixed viscosity ratio, we can easily show that there are interfacial positions, fa, for which A satisfies both problems. It is unclear whether in general this occurs and we have not analyzed much further. 144 Chapter 5 Static m u d channels As stated in Chapter 4,forcertain combinations of dimensionless parameters, it is possible for one of the fluids to be stationary on the narrow side of the annulus. Physically, this occurs when the modified pressure gradient is not large enough to overcome the yield stress of the fluid. From an industrial perspective the flows within this category that have most interest are those in which fluid 2 is static, either partly or fully, i.e. there is a static mud channel. Since the stability studies only predict linear stability of theflow,we ask whether a nonlinear perturbation could shift the mud channel. We study the displacement via numerical simulation. We will solve numerically the full system, V-[pV* ] t = - V - [ S + f], (5.1) where (5.2) |V*|=0 <S> |S|<7y/i?. tt((U,t) = o. *(1 ,t,t) = Q(t) (5.3) (5.4) (5.5) *(0,O,t) = * i „ f > , t ) , (5.6) *(4>,Z,t) =* (<t>,t), (5.7) mt (5.8) 145 Chapter 5. Static mud channels Whereas for the stability studies, we have Q(t) = 1, here we consider varyingflowrate, such as would •be possible with pressure pulsing. For the velocityfieldwe use an augmented Lagrangian algorithm, [26, 31], which fully describes the unyielded regions of theflow,the method is outlined in §5.1. We will work in afiniteelement approximation setting for the velocityfield,with bilinear basis functions and square elements, for an summary of thefiniteelement method see §5.3. The concentration equation is fully discretized using a conservative finite volume approximation at the middle point of each element and solved using a Flux Corrected Transport scheme, see §5.4. The FCT scheme is used to minimize both numerical dispersion and diffusion. Finally in §5.5 we investigate the effect of pulsing the flow rate on the displacement fronts. 5.1 Augmented Lagrangian Method The augmented Lagrangian method, although it accurately predicts unyielded regions in the flow of a visco-plastic fluid, is not widely used for simulating this kind of flow. Instead a regularised version of the constitutive laws of thefluidis used. This consists of considering, at the point of singularity, (i.e. shear rate equal zero),that the fluid has instead a very large Newtonian viscosity. For some recent work where the augmented Lagrangian algorithm has been successfully implemented we refer the reader to [80, 81, 82]. This section mainly gives an outline of the augmented Lagrangian algorithm, we will follow closely [26, 31]. 5.1.1 Principles of the method For the purpose of this thesis we will limit ourselves to afinitedimensional problem, i.e. since in any case we eventually discretise and solve afinitedimensional problem. Augmented Lagrangian As an example, let A be a symmetric, positive definite N x N matrix and suppose that b G R . We associate the following quadratic functional to these parameters J(y)= -{Av,v)-{b,v) l 146 (5.9) Chapter 5. Static mud channels where in (5.9), (•, •) is the canonical Euclidian inner product in R^. Let T be a linear mapping from R into R . Now, consider the minimisation problem N M J(u) < J(v) W e KerT u e (5.10) KerT. (5.11) Because of the properties of A it is clear that (5.10) has a unique solution, but the constraint (5.11) may make the solution difficult to compute. Introducing a Lagrange multiplier p to (5.10), in order to get an unconstrained problem, we have min v e R (5-12) {J(v) + (p,Tv)} N Now the problem is unconstrained but the Lagrange multiplier p appears as an extra unknown. This can be obtained through the solution of a saddle-point problem. Thus, we define (5.13) L(v,q) = J(v) + (q,Tv) recall that (u, p) will be a saddle-point of L if L(u,q)<L(u,p)<L(v,p) (5.14) Vv€R, ,qeR N M thus we have min max L(v, q) = max min L(v, q) veR qeR. N eR. veR M M = L(u, p) N q (5.15) It can be shown that L admits at least one saddle-point (u,p) on R x R , where u is the N solution of (5.10) and is common to all the We define the augmented Lr(v,q) Lagrangian L r M of L on R x R . saddle-points w M , for r > 0, by = J(v) + (q,Tv) + ^\Tv\ 2 147 = L(v,q) + ^\Tv\ 2 (5.16) Chapter 5. Static mud channels and note that a saddle-point of for which L(v,q) Ker u G T, will also be a L (v,q). However, (5.16) may be better conditioned to solve numerically. 5.1.2 ALG2 r saddle-point of In [26, 31] the authors describe an algorithmforfinding the saddle-points of an augmented Lagrangian functional L , they call it ALG2. r Let us consider the minimization problem mm{F(Tv) (5.17) + G(v)} v G V as our generalized minimization problem, where • V and H are Hilbert spaces. • T : V ^ H . • where F and G are convex, proper, and lower semi-continuous (l.s.c) on H and V, respectively. It is know that (5.17) and min {F(q) + G(v)} - (5.18) (v,q)€W with W = {(viq)€VxH,Tv-q = 0}. are equivalent problems. Introducing the variable q, linked to v through the linear equality equation Tv = q we can reformulate (5.18) as a saddle-point problem. We will in fact take the spaces V and H to be V = V p and H respective inner products and associated norms. We define, for v = L (il) 2 £ V fl, q € H Lagrangian functional L(v,q,p)=F(q) + G(v) + 148 x L (ft) with their (p,Tv-q) p 2 and p G H the Chapter 5. Static mud channels and then for r > 0, the augmented Lagrangian functional is defined as L{v,q,p) = L(v,q,p) r + (5.19) -\\Tv-q\\ . r 2 The iterative method A L G 2 that we use below is in fact a method for calculating the saddlepoints (u,p,X) of (5.19). ALG2 . (p^A ) €H xH arbitrarily specified; (5.20) 1 with (p , A") known, determine succesively n_1 G(v)-G(u ) + n +r(Bu u n F{q) Va n - (A , n n+1 by n - V)) > 0 Vv x H, and q-p ) + r(p n p n x G H n => Convergence of n l yn+l n 5.1.3 n G Vp, (5.21) G Vp - F{p ) G H n n (X ,T{v-u )) -p ~ ,T(v n u ,p ,X + Tu 11 n (5.22) H. ( n pn - Tu , q-p ))>0 _ ny p (5.23) ALG2 Results on convergence of A L G 2 are presented below, we just give a summary of what conditions T, F and G have to satisfy and state the convergence result. For a full description of the proof we refer to the reader to [26]. 149 Chapter 5. Static mud channels Conditions for T , F and G (i) T should be injective and the image of T should be closed in H. (ii) = 00 lim ^ \Q\->°° (5.24) \q\ (iii) Because G is convex, proper and l.s.c and using (i) and (ii) we have that (5.17) admits a unique solution u. (iv) Let us assume that F = FQ + F\, where Fi is convex, proper and l.s.c. on H, and Fo is convex, single valued, Gateaux differentiable on H, and uniformly convex over the bounded subsets of H in the following sense: For any M > 0, there exists a continuous function 6M'•[0, 2M] — > R, strictly increasing, with SQ = 0, such that for all p, q £ H, \p\ < M , \q\ < M , we have: (5.25) (Fo (q)-F (p),q-p)>5 (\q-p\). , , 0 M Now we are in a position to state the main result for convergence of the algorithm. Theorem 3 We assume that L assumption on T, F, G (i)-(iv) r admits a saddle-point and if p the following onV xHxH. p Under the satisfies n convergence results u —> u p —> p n n (5-26) (^Y^)r, 0<Pn=P< we have for ALG2 (u,p,X) 150 strongly strongly in V, (5.27) in H, (5.28) Chapter 5. Static mud channels n+l _ n _^0 X X A" i strong y i is bounded n H ( 29) i 5 (5.30) in H. Proof For a proof see [26]. • 5.2 5.2.1 Application of A L G 2 to the displacement problem. Steady model Let us take T(u) = Vu, G(u) = 0, F(q) - F (q) + F (q), 0 i where F (q) , /1 \VV+q\* , 1/2) = J UJ ^T-te Fi(q) = f 0 I^|W* + g|dft, \ + q-ndn, q € H x H, (5.31) qeHxH. Let us stop for a moment to make a couple of remarks. We have homogenized the boundary conditions with \I>* defined in (3.13). Recall that if ./HQO / Owe are in a Hilbert space setting, for the limit pm —> 0 we would have to work in the Banach spaces W ' * '™^) 1 1 1 for $ and (L ' / (J7) x L^ 1 1+1 m 151 / {9)) l+l m for q. The Chapter 5. Static mud channels augmented Lagrangian algorithm requires the solution space to be a Hilbert space. Nevertheless, once we have discretized our variational problem using afinite-dimensionalnumerical method (Finite Element Method for example), the approximate solution space V ^ is a finite dimensional p Hilbert space. The subscript h denotes the mesh scale and although the technical difficulty remains in the limit h —* 0, in practice we compute the approximate solution with a small finite h. Thus it is safe to assume that ourfinitedimensional solution space to be V h C Pt JY(f2) and 1 similarly we can take q G H. Outline of A L G 2 for a finite dimensional setting For a given constant r > 0 and using the definitions in (5.31) our augmented Lagrangian functional is, Lr(v,q,n)=F(q)+ f Jn • (Vv - fi J / |Vv - q)dft + Jn 1 Therefore, if (u,p, A), where u e Vb, and u is defined as * = 2 + u, is a saddle point of L , p Vr > 0, then (5.32) q\ dfl. r + u is a solution to (3.25). Now we are in a position to apply A L G 2 to our problem. Each step of the algorithm consists offindingu , p , X n n n+1 by solving sequentially: / A"-V(u-u") + (rVw"-p"- )-V(7j-u")dft > 0, 1 Jn Vv G Vb,p, u e Vo, n F(q)-F(p )+ n f \ n • (q- p ) + r(p n n Jn — Vw ) • n p (q - p )dfl n Vq e H, p n X+ = \ n 1 n + (Vu -p ), n P n n (5-33) > 0, (5.34) e.H p n > 0. (5.35) Minimizing equation (5.33) with respect to v is equivalent to solving the following Poisson 152 Chapter 5. Static mud channels equation, rAu = rV • p ~ n n l - V • A , (5.36) n where the right hand side is known from the previous step. The second step consists of minimizing equation ( 5 . 3 4 ) with respect to q in order to obtain p . This second step is the critical n one and shows the advantage of this method. For the case of visco-plastic fluids, in this step we will fully describe the yielded and unyielded regions. In minimizing ( 5 . 3 4 ) with respect to q in a discrete setting, we may locally minimize the integrand of L (u , q, A " ) with respect to n r q G Hh x Hh where Hh denotes the restriction of if to a particular discrete element, thus p = mi{-J ^L -i n o 7 d + s g - f + ^ | V * * + q| + - | q | 2 - ( A " + rV U ")-q}, (5.37) equivalently we have, 1 p n = /.|v**+g| q 2j l / 1/2^ 2 i{ ^^ds+^|V**+g|+-|V**+q| -(A"+rV(r+«")-/)-(V**+ )}. 2 g (5.38) Clearly, this expression will achieve its minimum when (V** + q) || ( A + rV(** + u ) - f). n n Therefore, letting + rV(** + u")-/), (V**+p") = 0(A" x = |A +rV(**+u")-/| (5.39) n (5.40) we have to find the minimizer of r (0x)2 1 M { 0 ) = Y(S / ) J ^ 7 ^ 2j o d r TV 2 s + + 2( o 0 X ?- 0 X - ^ 5 4 1 ) If x < ry/H then 9 = 0 minimizes M(9), thus the fluid will be unyielded. If x > ry/H we solve numerically for 9 the following equation, (9x) X / + ^ F H + r 0 x - x = 0, • (5.42)' which has a unique solution and therefore here the fluid is yielded in this iteration. Finally we update A using (5.35). 153 Chapter 5. Static mud channels 5.2.2 Transient model The derivation of the augmented Lagrangian algorithm for the transient case is basically the same as the one for the steady model. The main difference is with the augmented Lagrangian functional defined as, L (v \q \n) n+ = F(q ) n+ + n+1 r f pi • ( W * + 1 - q n + 1 ) d f 2 + J / |Vi> Jn n+1 - q" | df2 +1 2 * Jn (5.43) where the super index v means v(t ) where t n+1 n+l Thus the transient version of A jT \ n + l k • V(v n+1 - u ) n +l k =t + At. n+1 n will be, L G 2 + ((r +• £)Vu£ - -LVttf +1 - p t\) • n k Vv n+1 F(q ) n+1 - F(p£ ) +1 + f \+ n Jn l • (q n+1 - p ) n +1 k + r(p£ +1 V(«" - G q ) Vo, , u n+1 Vq" + 1 +1 e V, (5.44) pJJ )dfi > 0, k - 0 p +1 G H, pl KX\=K +P^K -Pl ), +1 x n +1 p • (q n+1 )dQ > ' 0, +1 +l +1 GH (5:45) Pk > o. (5.46) the subscript k is just the step of A L G 2 . We will update A n + 1 using (5.46) and the minimiza- tion problem (5.45) corresponds to finding the roots of (5.42). The only difference is in minimizing (5.44), which is equivalent to solve for each time step the following Poisson problem, (r + ^ ) A < + = r V - p - - V - A " - ^ V ^ . 1 n 1 (5.47) We solve (5.36) & (5-47) using a finite element method. In the following section we give a brief •outline of the method, which closely follows [43]. 154 Chapter 5. Static mud channels 5.3 Finite Element discretization In §5.2 we stated that the minimization problem (5.33) is equivalent to solve the Poisson problem (5.36), thus the variational problem for (5.36) is, Find u £ Vb, such that p a(u,v) = b(y) (5.48) V v £ Vo, p where a(u, v) = b(v) / Vu • Vudfi, • Jn = [ Jn fvdtl, with / being the right hand side of (5.36). The finite element discretization consists in replacing VQ by a finite dimensional subspace IP Vh, i.e. instead of solving (5.48), we solve Find u h £ V h such that a(u ,v) h = b(v) (5.49) Vv£V , h this approach is also known as the Galerkin method and the finite-dimensional subspace V/j is called an ansatz space. As an example, let us take to be a partition of our domain ft into closed triangles K with the following properties "1- V = U K e T h K , 2. For K, K' £T , K ^ K' and interior^) f| interior^') = 0. h 3. If K ^ K' but K f| K' = 0, then K f) K'is either a point or a common edge of K and K'. Let us define h as h = max{diam(/Y)|/Y e T }, h K is called an element of the partition and the vertices of the triangle are called the nodes. We call this a linear finite element approximation of our problem if 14 is defined as follows, V h = {u £ C(U)\ u\ K £ P\{K) . 155 \/K £ Th, u = 0 on dfl,} Chapter 5. Static mud channels where P\ (K) a + Bx + 7y is the set of polynomials of first degree (in 2 variables), i.e. pG P\{K) & p(x,y) = V (x,'y) G K and for a, 3, 7 G R. Now, let us consider the local interpolation problem p(oj) = iij for i = 1, 2, 3 (5.50) where p G Pi(K) and ai, a , 0 3 are the vertices of K, thus 2 p(x, y) = ui<pi(x, y) + u if2(x, y) + U3ip (x, 2 3 y), where tpi(aj) = 5ij, has a unique solution, see [43] and that U\K can be composed continuously, i.e. if px G P\{K) is the unique solution of (5.50) then u(x) = PK{X) for x G K G T. n Stability and convergence results for the finite element method for linear triangular elements can be found in [43]. Here we just state the result. If the weak solution u of (5.36) has weak derivatives of second order, then \\u - Uh\\i < Ch, where C depends on u but not on ft. We use this type of linear elements with linear basis functions ^ for solving numerically system (2.63) in §3.3.3. For the purposes of this chapter, we choose to work with quadrilateral elements and bilinear basis functions. All results stated above can be extended to this type of finite element discretization, (see [43] chap. 3 for reference). 5.4 Flux Corrected Transport Algorithm Now we turn our to attention to equation (5.8). Note that this concentration equation is coupled with system (5.1)-(5.7) via the stream function. Thus we will solve the full coupled 156 Chapter 5. Static mud channels system (5.1)-(5.8), using the augemented Lagrangian method in a finite element setting for the stream function and a Flux Corrected Transport (FCT) algorithm in afinitevolume setting for the concentration. We describe the FCT method below. For a full description of the algorithm we refer to the reader to [89]. Here we just give an outline of the method applied to our particular problem. It is well know that high order schemes present numerical dispersion near discontinuities and regions of high gradients, which translates into the growth of unphysical oscillations in the conserved quantity, (here c). In contrast, low order schemes do not suffer from this dispersion, but numerical diffusion is introduced in to the solution. Flux limiting methods, such as Flux-limiters and FCT have developed over the years to minimize both numerical dispersion and diffusion. We will concentrate on FCT schemes which limit thefluxby minimizing the overshoots and undershoots in the higher order numerical solution. The discretization will be made using a finite volume method which is well known to be conservative. The main reason that we use quadrilateral elements with bilinear basis functions for the stream function is that we will calculate the concentration field at the middle of each element, thus the calculation of numericalfluxesacross the boundary of each element is easier to obtain than with triangular elements. Let us write equation (5.8) in the following form, ldC df dg e dt dq> dt, , . where c = He, f = vHc = 9 = wHc = -> Calculating the numerical transportfluxesacross each element and using a forward Euler scheme 157 Chapter 5. Static mud channels for time we have: /~m+l At , „ /-m j-, n At, s n n . n (5.52) where F and G are the discrete numerical fluxes. Note that (5.52) is a second order approximation in space for C, thus numerical dispersion will be introduced in the solution, thus we use a FCT scheme to minimize this effect. First we calculate F L . and G L r 1 using a low order scheme in the following way, Vi n .CJ , . >0 HV , 1 N I V7! x .CS.!,- (5.53) i f V " ! . <0 and w.j.iCj ,- i f i y p . ^ j ^ o 1 G r . . c ^ i if w (5.54) ™ , <o. Now we compute the low order time advanced solution, A *-£ F , 1 .) — x 7 ( C f . , i — G . ?.— ± . 1 ' 0- L AC \ 7.1-1-4 *).? — (5.55) ~ Then we compute F^fi . and G . i using a high order scheme in the following way ff (5.56) (5.57) We compute the antidiffusive fluxes, to minimize the diffusion introduced by using a low order scheme, (5.58) tun K h i^- = +G (5.59) + We limit this antidiffusive fluxes to minimize the dispersion effect that the high order approximation has on the solution, (5.60) (5.61) A. 1 158 Chapter 5. Static mud channels and finally we apply the limited antidiffusivefluxesto the get the final solution, (5.62) -( The only thing that remains unclear is how the corrections L i j i+ and L .j i t + will be chosen to limit C^f so that it does not exceed a maxim value Cj"? nor a minimum value CJf . This 1 x n will be achieved in the following way. First we define the following quantities, 1. We compute the sum of all antidiffusivefluxesinto the grid point (i, j) (the centre of our element), P± = max(0,^_i j) - mm(0,A ij) max(0, + i+ ) - max(0,^. + i). 2. We compute the sum of all antidiffusive fluxes away from the grid point (i,j), FT. = max(0, A i + l .) - min(0, A^x-) + max(0, A i j + i) - max(0, A _i). id 3. Define s-rmax _ i,j — m a v m m (r<n rm V°i,j' rm rm rm \ °i-lj> 4. Define ( = qnox _ C5)A0AC, (C<° - C™")A0A£. 5. Define min(l,Q+/P+) if/ft >0 0 ifP+=0 min(l,Q-./F--)_ if P". > 0 0 159 (5.63) (5.64) Chapter 5. Static mud channels 6. Finally, the corrections become nan(Rt j,Ri ) +1 (5.65) d min(i?+.,i?- .) +lj i f V * J < min(i?+. ,i? ) if A-, min( R+.,i?- ) if A - , +1 J ij j+1 0 (5.66) Now, we are in the position to solve the system of equations (5.1)-(5.7) numerically and investigate the effects of pulsation of the flow rate in the displacement flow. 5.5 Numerical simulations for the annular displacement with pulsating flow rate In this section we are interested in the behaviour of the displacement flow in the presence of a pulsatingflowrate, thus (5.5) becomes (5.67) #(!,£, t) = 1 + dpsincot, Thus we will solve system (5.1)-(5.7) with (5.67) instead of (5.5), using the augmented Lagrangian algorithm in a finite element setting for the velocityfieldand for the concentration we will use the FCT scheme described in the previous section. 5.5.1 Validation To benchmark our code we follow the approach in [58]. For small eccentricities and certain combinations of the physical parameters, it is possible tofindan analytical solution that describes the shape of a steadily advancing displacement front, see [58]. In [58] it is shown that the c = 0.5 level set, viewed in a moving frame of reference, converges to the steady profile. Here we have performed similar benchmarking computations andfinalsimilar results. Figure 5.1 shows an example of these computations. We solve the steady model for the velocity, so that the only time dependence is through the concentration equation, i.e. 5 = 0 in (5.67). The P following physical and rheological parameters are used, K\ = 0.5, K = 0.4, m\ = 1, m = 1.2, 2 2 p = i p = 0.9, 7V,i = 1, 7V,i = 0.9, 7Y,! = 0.7 e = 0.0 and 3 = 0. Results are consistent x ) 2 160 Chapter 5. Static mud channels with the ones found in [58]. The validity of our code for the transient velocity field is shown in §3.3.3. Assuming fixed concentration we have seen the decay of the transient velocity to the steady state, at constant flow rate (5 = 0). P 0.8 0.6 0.4 0.2 N 0 -0.2 -0.4 -0.6 -0.8 0 t (a) t (b) (c) Figure 5.1: Convergence of computational solution to stable steady state displacement: a) Successive contours of c(</>, z, t) = 0.5 at times t = 0, 0.2, 0.4, ... 6. b) Analytical solution for the steady displacement, c) Convergence of the wide and narrow side interface positions to the steady state. Rheological and physical parameters: «i = 0.5, K = 0.4, m\ = 1, m = 1.2, =l = 0.9, 7y,i = 1, ry,i = 0.9, TY,I = 0.7 e = 0.0 and (3 = 0. 2 P l 2 ; P 2 5.5.2 Pulsation effects Pulsation at the beginning of the displacement In Figure 5.2 we show the initial condition for all the following displacements. One third of the eccentric Hele-Shaw cell is full of displacing fluid, cement, the remaining two thirds is full of the displaced fluid, mud. We will always begin with a flat interface, perpendicular to the £ axis between the fluids. For the rest of the section we will fix the rheological and some physical parameters to: K\ = 0.5, K 2 = 0.4, mi = 1, m 2 = 1.2, p\ = 1, p 2 = 0.9, Ty,i = 1, ry i t = 0.9, 7y,i = 0.7, (3 = 0. We will study the effects of eccentricity, the amplitude and phase of the pulsation and e on the displacement fronts.- 161 Chapter 5. Static mud channels •y 1.5 1 0.5- 0.2 0.5 0.8 Figure 5.2: Initial condition for the displacements. In Figure 5.3 we show the unsteady displacement of cement displacing mud. The eccentricity is e = 0.3. As we can see from the velocity field plots, both fluids are fully yielded along the annulus. We show times T/4, T/2, 3T/4 and T, where T is the period of the pulsation. For the steady state velocity model is defined as T defined as T = 2TT./U)€. = 2-K/U), for the transient velocity model is Recall that velocity changes happen on a fast time scale. Thus, in order to make a comparison of our models, we will work with the pulsation on the same time scale as the concentration. Figure 5.4 shows the effects of the pulsation in the steady model, we increase the amplitude to 5 = 0.1. No major effects can be seen on the displacement front and P velocity fields, as expected, both fluids are fully mobile along the annulus. In Figures 5.5 - 5.7 we investigate the effects of the pulsation in the transient model. Also we are interested in the results of varying e, recall that e is the ration between the advective and the viscous time scales. For Figure 5.5 we set e = 0.1 and u> = 10, note that if e —* 0, we recover the pseudo-steady model, this is consistent with our findings where we can see that the effects of the pulsation on the displacement front is minimal as in previous figures. If e ~ 0(1), viscous and advective time scales have the same order, thus the perturbation is picked up by the transient behaviour of the velocity as shown in Figure 5.6. Figures 5.7 - 5.9 show the reversed case, i.e. where viscous effects happen in a slow time scale, e = 10. In Figure 5.7 shows the case when the 162 Chapter 5. Static mud channels perturbation of the flow rate is fast, UJ = 10. Because we show the displacement front over the full period of the pulsation, time is not large enough for viscous effects to take place. When ui = 1, Figure 5.8, time is large enough to see viscous effects and the unsteady displacement takes place. Figure 5.9 shows the effects of a perturbation of theflowrate over a long time, we see clearly the viscous effects had taken place and a very long finger has developed along the annulus. Note that for Figures 5.8 & 5.9 we solved the problem in a moving frame of reference, the finger grows so long that we were unable to consider a domain large enough to capture the full finger. For Figures 5.10 - 5.12 we increase the eccentricity to 0.5. As shown in [58], as eccentricity increases we attain displacements that burrow slowly up the narrow side of the annlus, i.e. although the interface is propagating along the annulus, in the far field the velocity field is unyielded in the narrow side. In Figure 5.10 we got the same results as the ones found in [58], this was expected as we solve the steady-state velocity model. Figure 5.11 shows the effects pulsation has on the steady state model. We fix the amplitude of the pulsation, S = 0.1. As p above, we cannot see major changes in the displacement front but an interesting phenomena is the decrease of the unyielded regions on the far field, as is clearly seen in the velocity field plots. On the other hand, Figure 5.12 shows that the transient model actually captures the effects of the pulsation. When we reach the maximum point on the pulsation, Figures 5.12 a)-b) bothfluidsare fully yielded, even in the far field, as we decrease the intensity of the flow rate stagnant zones appear until we reach the minimum point of the pulsation and both fluids are unyielded on the narrow part of the annulus for the far field. When the intensity of the flow increases again we fully remove the stagnant zones. One interesting phenomena is that in contrast to the steady state model, the velocity field of the transient model does not show the same burrowing motion. In Figures 5.10 and 5.11 it is clear that in the vicinity of the interface the displacingfluidis forced to the narrow side and the displacedfluidis forced away from the narrow side. For the transient model this burrowing motion is expanded along the annulus- not only near the interface. This will lead us to investigate the effects of the pulsation technique and transient effects of the velocity field in the removal of a mud channel. Figure 163 Chapter 5. Static mud channels 5.13 shows the predicted mud channel width using the steady model without pulsing the flow rate. Figures 5.14 and 5.14 show the width of the mud channel with pulsation amplitudes 8 = 0.1 and 5 = 0.2 respectively. As in the previous cases there is no clear evidence of a P P change in the displacement. Figures 5.16 and 5.17 show the effects of a pulsating flow rate on the transient model with amplitudes, S = 0.1 and 5 = 0.2 respectively. Clearly the effects of p P the pulsation are fully captured by the transient velocity field, as above we have a burrowing motion expanding along the annlus, i.e. even though we have a mud channel in the narrow side of the annulus the far field velocity is clearly yielded outside this region. As we increase the magnitude of the pulsation, this motion expands further to the narrow side of the annulus, and therefore a decrease of the width of the mud channel is achieved. 164 Chapter 5. Static mud channels (b) (a) 4> (e) (d) (c) $ <f ' • (g) (f) (h) Figure 5.3: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2TX/U), U = 10 and S = 0. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, m\ = 1, m = 1.2, p\ = 1, p = 0.9, Ty,i = 1, 7y,i = 0.9, Ty,i = 0.7 e = 0.3, /3 = 0. p 2 2 165 2 Chapter 5. Static mud channels (b) (a) (d) (c) in HI ii'" flMI 0.2 (e) 0.5 0.8 (h) (f) Figure 5.4: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2n/uj, OJ = 10 and S = 0.1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p — 0.9, Ty,\ — 1, ry,! = 0.9, 7 y i = 0.7 e = 0.3, /3 = 0. p 2 2 166 2 Chapter 5. Static mud channels (a) (e) (b) (c) (f) (d) (g) (h) Figure 5.5: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 27T/W, UJ = 10, 5 = 0.1 and e = 0.1. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: KI = 0.5, « 2 = 0.4, mi = 1, m = 1.2, p\ = I, p — 0.9, Ty,i = 1, TY,I = 0.9, 7y,i = 0.7 e = 0.3, /3 = 0. V 2 167 2 Chapter 5. Static mud channels (a) (e) (b) . (d) (c) (f) (h) (g) Figure 5.6: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2-n/uie, u> = 10, 5 = 0.1 and e = 1. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, ry^ = 1, ry,i = 0.9, Ty,i = 0.7 e = 0.3, 8 = 0. P 2 2 168 2 Chapter 5. Static mud channels 0.2 (a) (e) (b) 0.5 0.8 0.2 (d) (c) (f) 0.5 (h) (g) Figure 5.7: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2n/uje, UJ = 10, S = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, Ty \ = 1, TY,I = 0.9, TY,I = 0.7 e = 0.3, /3 = 0. p 2 2 169 2 t Chapter 5. Static mud channels m 3.5 3 2.5 2 1.5 1 0.5 0.2 0.2 0.5 (a) 0.5 0.8 (d) (c) (b) 4HH r 2.5! - \ 0.2 (e) 0.5 (f) (h) (g) Figure 5.8: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2n/u>e, OJ = 1, S = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K% = 0.4, m\ = 1, m = 1.2, p\ = 1, p = 0.9, Ty,i = 1, TYI = 0.9, TYi = 0.7 e = 0.3, 8 = 0. p 2 170 2 Chapter 5. Static mud channels 0.2 0.5 0.8 4> . 0.2 0.5 0.8 <f (a) (b) Figure 5.9: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 27r/we, UJ = 0.1, 5 = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. Physical and rheological parameters: K\ = 0.5, = 0.4, mi = 1, m = 1.2, = 1, p = 0.9, Ty = 1, TY,I = 0.9, 7y,i = 0.7 e = 0.3, /5 = 0. p K 2 2 P l 2 A 171 Chapter 5. Static mud channels (a) (b) < > ) - • (d) «t> < > t <> (f) (e) (c) (h) (g) Figure 5.10: Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2-ir/uj, UJ = 10 and 5 = 0. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, Ty,i = 1, TY,I = 0.9, ry = 0.7 e = 0.5, (3 = 0. P 2 2 A 172 2 Chapter 5. Static mud channels ii J iiii n, I.,. m : lit III 'In In In In iii!: •ill: 0.5 • (a) 0.2 0.8 (b) 0.5 0.8 (d) (c) !!: in in " i n ii. -iii in in In III 1 III,, ii,;; •i! Ill III III 0.2 0.5 (e) 0.8 (h) (f) Figure 5.11: Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2TT/UJ, U> = 10 and 5 = 0.1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, K = 0.4, mi = 1, m — 1.2, p\ = I, p = 0.9, ry,\ = 1, Ty,i = 0.9, 7 y i = 0.7 e = 0.5, 8 = 0. p 2 2 173 2 Chapter 5. Static mud channels o o (a) o (b) (e) (c) <f 0 o «. (d) 0 « (g) (f) (h) Figure 5.12: Displacementflowin eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2ir/ue, OJ = 10, 5 = 0.1 and e = 0.6. Interface position and velocityfieldfor times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: KI = 0.5, «2 = 0.4, mi = 1,TO = 1.2, p\ = 1, p = 0.9, ry,i = 1, TY,I = 0.9, ry,! = 0.7 e = 0.5, (3 = 0. p 2 174 2 Chapter 5. Static mud channels Figure 5.13: Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Steady state velocity model. Period T = 2TT/U, to = 10 and S = 0. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K.\ = 0.5, K2 = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, Ty,i = 1, ry,i = 0.9, TY,I = 0.7 e = 0.8, /3 = 0. p 2 175 2 Chapter 5. Static mud channels ,-»•. J 41 0.2 (a) 0.5 0.2 (b) 0.5 0.8 (d) (c) ,,,, ""•Mil i, i i,i 1 in in In in in in 05 • n iii 0.2 0.5 0.8 0.2 (f) 0.5 (h) Figure 5.14: Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Steady state velocity model. Period T = 2n/u, u = 10 and 6 = 0.1. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: «i = 0.5, K = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, Ty,i = 1, Ty,i = 0.9, 7 v , i = 0.7 e = 0.8, /3 = 0. P 2 2 176 2 Chapter 5. Static mud channels in' 2.5:; 2.5fc ... ^ 1.5 ^ 1.5 *M.5 0.5h 0.5 [ *1 0.2 (a) 0.2 0.5 0.8 (b) 0.5 0.8 0.2 0.5 0.8 (d) (c) in* 2.5 2.5;, I 1.5 "J-1.5 •1.5 In In In 0.5 0.5 j! 0' (e) 0.2 0.5 0.8 0.2 0.5 0.8 (f) (g) (h) Figure 5.15: Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Steady state velocity model. Period T = 2-rr/uj, LO = 10 and 5 = 0.2. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, KI = 0.4, m\ = 1, mi = 1.2, p\ = 1, pi = 0.9, Ty,i = 1, Ty,i = 0.9, ry,! = 0.7 e = 0.8, /3 = 0. p 177 Chapter 5. Static mud channels 8811 0.2 0.5 0.8 0.2 (a) 0.2 (b) 0.5 0.5 0.2 (c) 0.8 0.2 (e) 0.8 0.8 (d) 0.5 0.8 (g) (f) 0.5 (h) Figure 5.16: Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Transient velocity model. Period T = 2ir/coe, u> = 10, 5 = 0.1 and e = 0.6. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K,\ = 0.5, K = 0.4, mi = 1, m = 1.2, pi = 1, p = 0.9, TY,I = 1, 7y,i = 0.9, 7 y , i =0.7 e = 0.8, 8 = 0. P 2 2 178 2 Chapter 5. Static mud channels (a) (b) • (e) (c) (d) • <> t (f) • (g) (h) Figure 5.17: Displacement flow in eccentric annulus, interface propagation. Mud channel formation. Transient velocity model. Period T = 2-jr/coe, co = 10, S = 0.2 and e = 0.6. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: K\ = 0.5, KI = 0.4, m\ = 1, rri2 = 1-2, pi = 1, P2 = 0.9, 7v,i = 1, TY,I = 0.9, 7Y,i = 0.7 e = 0.8, /3 = 0. p 179 Chapter 5. Static mud channels Pulsation after mud channel forms Now we turn our attention to the worst case scenario, when a mud channel forms in the narrow part of the annulus. We tried to simulate the case in which the mud channel has already formed and the flow may be described as a parallel flow in a Hele-Shaw cell with the fluid on the narrow side static. We use the same rheological and physical parameter as the previous section, here the interface between the fluids is placed at <f> = 0.8. Figure 5.18 shows the effects of the pulsation on the steady state velocity model with an amplitude of S = 0.2. This model predicts that p the mud channel will remain static. In Figure 5.19 we show the effects the pulsation has on the transient model with an amplitude of 5 = 0.2. Clearly as we go over the full period of P the pulsation the mud channel slowly begins to yield until it is fully moving, then goes back to the static mud channel after the pulsation period is over. Thus in each pulsation the mud channel will yield, move up the narrow side for a short period of.time and stop again. Another interesting behaviour is that it seems that the flow is stable, no growing instabilities were to be seen at the interface. Actually, the interface remained unperturbed for all time and thus the removal of the mud channel was unsuccessful. One of the reasons we can think of has to be with the choice of the concentration model to simulate this problem. Perhaps an interface tracking model is a better choice for this type of problems. Actually, we tried to solve the interface tracking model using a conformal mapping technique. Due to our transformation, we lost control over the hyperbolic problem defining the interface, allowing us to use only a low order method to solve the problem. We think that the numerical diffusion was to large and no effects of the perturbed flow rate could be seen. An adaptive technique could be a better approach to attack this problem. We cannot therefore be completely sure if pulsation of the flow rate, after the mud channel has already formed, helps its removal. Until we extend our approach to an adaptive technique we cannot derive any firm conclusion. This remains as future work and will not be investigated in this thesis. As a conclusion to the chapter, the transient model (5.1)-(5.7) with (5.67) instead of (5.5) fully captures the perturbations of the velocity field. This may lead to a reduction of the mud channel telling us that pulsation of the flow rate at the beginning of the displacement may be 180 Chapter 5. Static mud channels used as a tool to prevent or remove the mud channels along the annulus. •BMMHR * SBWHlf SillllllliilikE •••((••In •HI •SH 0.2 (a) 0.5 ^^^^^ 0.8 (c) (b) (d) (0 (e) Figure 5.18: Displacement flow in eccentric annulus, interface propagation. Mud channel. Steady state velocity model. Period T = 2-ir/cj, ui = 10 and 5 = 0.2. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: n\ = 0.5, K = 0.4, mi = 1, m = 1.2, pi = 1, p = 0.9, Ty,i — 1, r^i = 0.9, = 0.7 e = 0.8, 8 = 0. P 2 2 181 2 Chapter 5. Static mud channels (a) (e) ' (b) (c) (f) (g) (d) (h) Figure 5.19: Displacement flow in eccentric annulus, interface propagation. Mud channel. Transient velocity model. Period T = 2jr/ue, u = 10, 5 = 0.2 and e = 0.6. Interface position and velocity field for times: a)-b) T/4. c)-d) T/2. e)-f) 3T/4. g)-h) T. Physical and rheological parameters: « i = 0.5, n = 0.4, mi = 1, m = 1.2, p\ = 1, p = 0.9, ry,\ = 1, 7y,i = 0.9, P 2 2 Ty i = 0.7 e = 0.8, 3 = 0. t 182 2 Chapter 6 Summary and conclusions This thesis is concerned with modelling and analysis of various transient effects in primary cementing of an oil well. This process involves the displacement of one shear-thinning yield stress fluid by another along a narrow eccentric annular duct. 6.1 Main contributions of the thesis The contributions of this thesis fall into 4 categories: • We developed a transient version of the model found in [10] and demonstrated its validity for small e. • The nonlinear evolution equation for the velocity field (stream function) is a third order partial differential equation, in 2 spatial dimensions and in time. We proved that there exists a unique solution to this equation and that the solution changes continuously with the data, i.e. well-posedness. These results are detailed and are summarised in chapter 3. • In the case that a long finger develops along the annulus during displacement, we were able to analyse the linear stability of the limiting parallel multi-layer flows. Specific results are summarised below in §6.2. • In the case that a mud channel develops in the narrow side of the annulus, we demonstrated that by perturbing the flow rate, via pulsation, the width of the mud channel can be decreased. It is necessary to pulse the flow rate during displacement, rather than after the mud channel has already formed. Further results are summarised below in §6.3 183 Chapter 6. Summary and conclusions 6.2 Interfacial instabilities When a longfingerdevelops wefindessentially a multi-layer parallelflowalong the length of the annulus. In chapter 4 we studied the linear stability of theseflowsfor the case when bothfluidsare yielded at the interface. This stability problem depends on 14 dimensionless parameters in its full generality and hence very difficult to study analytically. Thus, we first consider two simpler cases, parallelflowof two Newtonian fluids and parallelflowof two power law fluids, both in a vertical Hele-Shaw cell with constant height, i.e. a concentric annulus. For both cases we found that buoyancy plays a major role in causing instability. In the absence of density difference it seems that theflowis linearly stable. Inclination of the Hele-Shaw cell has destabilizing effects, if the heavierfluidis on top of the lighter,fluid,theflowmay be unstable. Again, in the absence of buoyant forces, viscous effects and inclination are not strong enough to drive an instability. The stability problem for yield stressfluidshas to be solved via computation. Here we solve the full problem, i.e. the parallelflowof two yield stressfluidsin a vertical or near vertical eccentric Hele-Shaw cell. We consider the case when bothfluidsare yielded at the interface. If the eccentricity is equal to zero, i.e. a concentric annulus, wefindanalogous results to the Newtonian and power lawfluidcases: buoyancy is the key parameter for instability. Although buoyancy forces are the primary reason for instability to arise, viscosity differences play a secondary role in the instability. For example, for certain values of density and viscosity if Pi > pi and K\ > K2 theflowmay be stable, but if p\ > pi and n\ < K theflowmay be 2 linearly unstable. This shows that even though density differences drive the unstable behaviour, viscosity differences have stabilizing or destabilizing effects. This phenomenon is in contrast with the work of Gondret & Rabaud, [33]. They showed that in a horizontal Hele-Shaw cell, viscosity differences only determine the wave velocity of the growing instability. An interesting phenomena that is worth to mention is the inclusion of eccentricity in to the stability problem. Eccentricity can have stabilizing or destabilizing effects. For example, in the case of a concentric annulus we show that for a certain range of values of density and viscosity 184 Chapter 6. Summary and conclusions if pi < p2 and'/ti > K , the flow was linearly unstable for all wave lengths. If we increase 2 the eccentricity, the flow will be linearly stable for long wave lengths. In contrast, if the more viscous fluid is also heavier, the flow was linearly stable for all wave lengths for the concentric annulus. If we now increase the eccentricity theflowwill become linearly unstable for long wave lengths. From such results we inferred that all significant changes seem to happen for long wave lengths. Thus we considered the long wave length asymptotic behaviour of our problem. We presented several stability maps in §4.7. As expected for the concentric annulus, in the absence of buoyancy theflowwill be linearly stable. When eccentricity is different from zero, all the rheological parameters as well as buoyancy affect the instability regions. Thus, even without the presence of buoyancy forces and inertia, the flow may be unstable. We do not believe that these instabilities have been identified before. In [59] the authors use a lubrication model to predict whether a good displacement takes place or not. Surprisingly, our unstable domain occurs only as a subset of the unsteady displacement domain predicted by [59]. This means that if a long finger develops during the displacement, for certain parameter ranges the parallelflowmay be unstable. Therefore, in unsteady displacements this suggests that mixing may occur at the interface of the finger due to the instability that we study here. It is not clear why there is apparently no interfacial instability for parameter regimes for which the annular displacement would be steady. Note that for such parameters, a base parallelflowis a perfectly well-defined solution of the two fluid problem. 6.3 Mud channel removal In chapter 5 we studied the effects of a pulsatileflowrate for the removal of a mud channel. This is via numerical simulation. In chapter 4 we showed that if one of thefluidsis unyielded at the interface theflowwill be linearly stable. We investigated the effects that pulsation has on the mud channel byfirstusing the steady state velocity model. If we begin the pulsation after the mud channel already formed, this model predicts that the mud will remain static. On the other hand, if we use our transient model, 185 Chapter 6. Summary and conclusions it predicts that as we go over one period of the pulsation, the mud channel slowly begins to yield until the mud is fully mobile and goes back to be static after the period is over. Thus, in each pulsation, the mud channel will yield, moves up the narrow side of the annulus for a short period of time and then stop again. If instead we pulsate the flow rate from the beginning of the displacement, the transient model predicts that a reduction of the mud channel width takes place. The steady state velocity model with a pulsatile flow rate shows no mayor changes with respect to the steady state velocity model with constant flow rate. Therefore, the transient model is necessary to fully capture the effects of the perturbation of theflowrate. 186 Bibliography F. ANDREU, C. BALLESTER, V. CASELLES, J.M. Total Variation a Solute in Pulsating J. AZAIEZ A N D B. SINGH, Stability Proc. Roy. Soc. Flow Trhough a Tube, Proc. Roy. Soc. of miscible displacements of shear thinning interfaces N . RASHIDNIA, T. in a cylindrical tube . Rocks. A. BARNES, J.F. HUTTON (1989). G.K. BATCHELOR, (1989). Drilling BEIRUTE Muds in a MAXWORTHY A N D J. K U A N G , Instability of of Fluid Flows Through to Rheology. Elsevier Kluwer Academic Publishers (1990). H. R.M. fluids Phys. Fluids 17, 052103, pp. 1-11, (2005). G. I. BARENBLATT, V.M. ENTOV A N D V.M. RYZHIK, Theory Natural the Phys. Fluids. (14)5, pp. 1557-1571, (2002). cell. BALASUBRAMANIAM, miscible for through a Tube, of a Solute in a Fluid Flowing R. ARIS, On the Dispersion of A (259), pp. 370-376, (1960). R. Problem Journal Functional Analysis,180, pp. 347-403, (2001). Flow, R. ARIS, On the Dispersion A (235), pp. 67-77, (1956). Hele-Shaw M A Z O N , The Dirichlet AND An Introduction A N D R.W. by Cement K. WALTER, to Fluid An introduction FLUMERFELT, Mechanics Slurries Cambridge University Press Dynamics. of Using an Accurate the Displacement Rheological Model, Process of Society of Pe- troleum Engineers, paper number SPE 6801, (1977). S.H. B i T T L E S T O N , J. F E R G U S O N ment during primary an eccentric J.W. annular B R I C E AND Flow Techniques, cementing Hele-Shaw AND I.A. F R I G A A R D , Mud removal and cement place- of an oil well; laminar cell. R . C . H O L M E S , Engineered Casing Cementing of Multiple Casing, C. R. C L A R K A N D L . G . CARTER, Mud Displacement Petroleum Engineers, paper number SPE 4090, (1973). COUTURIER, displacements in Programs Using Turbulent J. Pet. Tech., pp. 503-508 (1964). M.A. CHILDERS, Primary Cementing paper number SPE 1919, (1968). M. non-Newtonian J. Engng. Math., 43, pp. 229-253, (2002). D . GUILLOT, H . HENDRIKS A N D F. ated Spacer Properties for Optimal Mud Removal Society of Petroleum Engineers, with Cement Slurries, C A L L E T , Design in Eccentric Annuli, Rules Society of and Associ- Society of Petroleum Engineers, paper number SPE 21594 (1990). P. COUSSOT, Mudflow Rheology and Dynamics, P. COUSSOT, Saffman-Taylor 363-376, (1999). instability A.A. Balkema Rotterdam (1997). in yield-stress 187 fluids. J. Fluid Mech. (380), pp. Bibliography [17] B. D A C O R O G N A , (1989). [18] P.G. [19] D. Direct Methods DRAZIN, Introduction DUSTREHOFT, Pulsation in the. Calculus to Hydrodynamic Cambridge Univeristy Press (2002). Stability. G . WILSON A N D K . NEWMAN, to Control AMS Vol. Springer-Verlag of Variations. Field Study on the Use of Cement Society of Petroleum Engineers paper number SPE Gas Migration. 75689 (2002). [20] G . DUVAT A N D [21] M.J. J.L. LIONS, Inequalities in mechanics I. E K E L A N D , R . T E M A M , Convex [23] P. E R N , F. CHARRU A N D P. LUCHINI, Stability stratified EVANS, Partial P. [27] Analysis Differential analysis of a shear flow with fluids strongly PETITJEANS A N D E. MEIBURG, Desntiy-driven in a Hele-Shaw cell. unsta- J. Fluid Mech. 451, pp. 239-260, (2002). Lagrangian North-Holland, (1983). Methods. LA. F R I G A A R D , S.D. HowiSON A N D I.J. SOBEY, On the Stability Fluid. SIAM (1999). Problems. American Mathematical Society, (1998). FORTIN A N D R . GLOWINSKI, Augmented a Bingham of Poiseuille Flow of J. Fluid Mech. 263, pp. 133-150, (1994). [28] I.A. FRIGAARD, Super-stable parallel Fluid Mech., 100, pp 49-76, (2001). [29] and Variational Equations. KUROWSKI, P. ble flows of miscible [26] M. B. J. Fluid Mech. 496, pp.295-312, (2003). viscosity. J. FERNANDEZ, IN: 'E. Schlumberger Educational Services (1990) Chapter 5. Well Cementing, [22] [25] Springer-Verlag (1976). ECONOMIDES, IMPLICATIONS O P CEMENTING O N W E L L PERFORMANCE. NELSON ( E D . ) , [24] L.C. and physics. flows D . GILBARG, N. S. TRUDINGER, Elliptic of multiple partial viscoplastic differential fluid. equations J. Non-Newt. of second order. Berlin ; New York : Springer-Verlag, (1983). [30] R. GLOWINSKI, J.L. inequalities. [31] [32] LIONS A N D R . TREMOLIERES, R . GLOWINSKI, Augmented Lagrangian of boundary-value North-Holland (1983). P. Numerical analysis of variational North-Holland (1981). problems. GONDRET, N. Methods: RAKOTOMALALA, M. Applications RABAUD, D . parallel flows in finite aspect ratio Hele-Shaw to the numerical SALIN A N D P. cell: Analytical solution WATZKY, and numerical Viscous results. Phys. Fluids 9(6) (1997). [33] [34] P. GONDRET A N D M. RABAUD, Shear,instabillity cell. Phys. Fluids 9(11) (1997). V.A. GORODTSOV A N D V.M. Newtoinan [35] fluids in a Hele-Shaw R . GOVINDARAJAN, Effect YENTOV, cell. of two-fluid Instability of the parallel flow in a displacement fronts Helle-Shaw of non- J. Appl. Maths. Mechs. (61)1, pp. 111-126 (1997). . of miscibility on the linear instability Int. J. Multiphase Flow. (30), pp. 1177-1192, (2004). 188 of two-fluid channel flow. Bibliography [36] N. G O Y A L A N D E. M E I B U R G , Unstable density stratification Hele-Shaw cell: influence of variable viscosity of miscible fluids in a vertical on the linear stability. J. Fluid Mech. 516, pp. 211-238, (2004). [37] Rheology J . HARRIS, and Non-Newtonian flow. [38] E . J . H I N C H , A note on the mechanism shearing [39] fluids. J. of the instability E . J . H I N C H A N D F. P L O U R A B O U E , Kelvin-Helmholtz G . M . HOMSY, at the interface between two Fluid Mech., 114, pp. 463-465, (1984). effect from the small region near the meniscus. [40] Longman (1977). Viscous Fingering instability in a Hele-Shaw cell: Large Phys. Fluids. (17)052107, pp. 1-13, (2005). Ann. Rev. Fluid Mech. (19), pp. 271- in porous media. 311, (1987). [41] A.P. H O O P E R A N D W . G . B O Y D , Shear flow instability fluids. [42] at the interface between two viscous Fluid Mech., 179, pp. 201-225, (1987). J. D.D. J O S E P H A N D Y.Y. R E N A R D Y , Fundamentals of Two-Fluid Dynamics, Part I & II. Springer-Verlag New York, Inc. (1993). [43] P. K N A B N E R A N D L. A N G E R M A N N , Numerical Differential [44] L. K O N D I C , P. P A L F F Y - M U H O R A Y A N D M . J . S H E L L E Y , Models Shaw flow. [45] Partial of non-Newtonian Hele- L. K O N D I C , M . J . S H E L L E Y A N D P. P A L F F Y - M U H O R A Y , Non-Newtonian instability. E. L A J E U N E S S E , in a Hele-Shaw J . MARTIN, cible displacement N. in a Hele-Shaw Hele-Shaw flow Physical Rev. Lett. (80)7, pp. 1433-1436, (1998). E. L A J E U N E S S E , J . M A R T I N , N. R A K O T O M A L A L A A N D D. cilbe Displacements [47] and Parabolic Physical Rev. E. (54)5, pp. R4536-R4539, (1996). and the Saffaman-Taylor [46] Methods for Elliptic Springer-Verlag New York (2003). Equations. Cell. S A L I N , 3D Instability of Mis- Phys. Rev. Lett. (79)26, pp. 5254-5257, (1997). R A K O T O M A L A L A , D. eel at high rates. J . S A L I N A N D Y.C. YORTSOS, Mis- Fluid Mech. (398), pp. 299-319, (1999). [48] A. L I N D N E R , P. C O U S S O T A N D D. B O N N , Rev. Lett. (85)2, pp. 314-317, (2000). Viscous [49] H.P. L A N G T A N G E N , Heidelberg (1999). Differential [50] CF. Predict [52] Partial C F . L O C K Y E A R A N D A.P. H I B B E R T , Integrated Factors for Field [51] Computational R.H. Succes, J . LOCKYEAR, D.F. and Prevent, fingering Equations Primary Springer-Verlag Berlin Cementing J. Study Defines Key Pet. Tech. pp. 1320-1325, (1984). R Y A N A N D M.M. G U N N I N G H A M , Cement Channeling: How to Society of Petroleum Engineers, paper number SPE 19865, (1990). M C L E A N , C . W . M A N R Y A N D W . W . W H I T A K E R , Displacement mary Cementing, Phys. in a yield stress fluid. Pet. Tech. pp. 251-260, (1967). 189 Mechanics in Pri- Bibliography [53] H . P A S C A L , Rheological interface [54] H . P A S C A L , Dynamics effect of non-Newtonian fluids on dynamics of moving Int. J. Engng. Sci. (22), pp. 227-241 (1984). of moving interface in porous media for power law fluids with a Int. J. Engng. Sci. (22), pp. 577-590 (1984). yield stress. [55] behaviour in porous media. H . P A S C A L , A theoretical for Bingham displacing analysis fluids of stability of a moving and its application interface in oil displacement in a porous medium mechanism. Can. J. Chem. Engng. (64), pp. 375-379 (1986). [56] C.W. P A R K A N D G.M. H O M S Y , The Fluids. (28)6, pp. 1583-1585 (1985). [57] [59] cementing displacements. cementing annuli: prediction computational simulation of eccentric IMA J. Appl. Math., 69, pp. 557-583, 2004. S. P E L I P E N K O & I.A. F R I G A A R D , Visco-plastic row eccentric in primary J. Engng. Math., 46(1), pp.1-26, (2004). S. P E L I P E N K O & I. A. F R I G A A R D , Two-dimensional annular Phys. of long fingers in Hele-Shaw flows S. P E L I P E N K O A N D LA. F R I G A A R D , On steady state displacements of an oil well. [58] instability of travelling fluid displacements wave solutions in near-vertical and interfacial nar- instability. J. Fluid Mech, 520, pp. 343-377, (2004). [60] F. P L O U R A B O U E A N D E.J. H I N C H , Fluids. (14)3, pp. 922-929, (2002). [61] R. R A G H A V A N A N D S.S. . medium. [62] Kelvin-Helmholtz instability M A R S D E N , The stability in a Hele-Shaw of immiscible cell. Phys. liquid layers in a porous J. Fluid Mech. (48), pp. 143-159, (1971). R. R A G H A V A N A N D S.S. flow of immiscible M A R S D E N , A theoretical liquids in a porous medium. study of the instability in the ^parallel Quart. Journ Mech and Applied Math. (26), pp. 205-216, (1973). [63] J.E. R I T T E R , A N D G.W. M C D A N I E L , Viscous Fluid Mud Displacement of Petroleum Engineers, paper number SPE 1772, (1967). [64] P.G. S A F F M A N A N D G. T A Y L O R , Hele-Shaw cell containing The penetration a mor viscous liquid, Techinque, of a fluid into a porous Society medium or Proc. Roy. Soc. A (245), pp. 312-329, (1958). [65] C.W. S A U E R , Mud Displacement During the Cementing Operation: A State of the Art, Society of Petroleum Engineers, paper number SPE 14197, (1985). [66] P. J. S C H M I D , D . S. H E N N I N G S O N , Verlag, New York, Inc., (2001). [67] and Transition J. S C O F F O N I , E. L A J E U N E S S E A N D G.M. H O M S Y , Interface ments of two miscible [68] Stability fluids in a vertical pipe . M. S H A R I A T I A N D Y.C. Y O R T S O S , porous media. Stability instabilities during Springerdisplace- Phys. Fluids. (13)3, pp. 553-556, (2001). of miscible Phys. Fluids. (13)8, pp. 2245-2257 (2001). 190 in Shear Flows. displacements across stratified Bibliography [69] R.C. S I M I T H , (1984). Succesful [70] T.R. S M I T H , 629, (1984). Cementing [71] J.P. STOKES, Primary Cementing Displacement D . A . WEITZ, J.P. Can Be a Reality, J . Practices-Field GOLLUB, C H A I K I N A N D H . M . L I N D S A Y , Interfacial in Axial Annular Flow, Int. J. M.O. ROBBINS, displacement P.M. in a porous Fluids in Eccentric Annular Geome- of One Newtonian Fluid by Another: Denstiy Multiphase Flow. (23)1, pp. 113-129, (1996). C.T. T A N A N D G . M . H O M S Y , Stability of miscible displacements in porous media: Recti- Phys. Fluids. (29)11, pp. 3549-3556 (1986). C.T. T A N A N D G . M . H O M S Y , Simulation of nonlinear viscous fingering in miscible dis- Phys. Fluids. (31)6, pp. 1330-1338 (1988). placement. [76] of immiscible Pet. Tech. pp. 564- J. Non-Newt. Fluid Mech, 45, pp 149-169, (1992). linear flow. [75] stability P. S Z A B O A N D O. H A S S A G E R , Displacement Effects [74] A . DOUGHERTY, P. S Z A B O A N D O. H A S S A G E R , Flow of Viscoplastic tries. J . [73] Applications, Physical Rev. Lett. (57)14, pp. 1718-1721 (1986). medium [72] Pet. Tech. pp. 1851-1858, G . T A Y L O R , Dispersion of Soluble Matter in Solvent Flowing Proc. Slowly through a Tube, Roy. Soc. A (219), pp. 186-203, (1953). [77] G . T A Y L O R , Conditions be Used to Measure [78] A . TEHRANI, placement under Molecular Which Dispersion Diffusion, of a Solute in a Stream S . H . B I T T L E S T O N & P . J . L O N G , Flow instabilities of one non-Newtonian fluid of Solvent can Proc. Roy. Soc. A (225), pp. 473-477, (1954). by another. during annular dis- Experiments in Fluis 14, pp. 246-256, (1993). [79] B . S . T I L L E Y A N D S . H . D A V I S , Linear channel. [80] strategy and some benchmark D . V O L A , F. B A B I K currents [82] G . VINAY, A . WACHS fluids. J. A N D J.F. waxy crude oil flows. J . unsteady flows of Bingham annulus. J. fluis: a Comp. Phys. (187), pp. 441-456, (2003). On a numerical strategy to compute gravity Comp. Phys. (201), pp. 397-420, (2004). AGASSANT, Numercial simulation of non-isothermal Non-Newt. Fluid Mech, 128, pp. 144-162, (2005). I.C. W A L T O N A N D S . H . B I T T L E S T O N , The axial flow of a Bingham eccentric plastic in a narrow Fluid Mech. 222, pp. 39-60, (1991). [84] S.D.R. W I L S O N , The Taylor-Saffman (220), pp. 413 425, (1990). [85] results. J . AND J.C. LATCHE, of non-Newtonian viscoplastic [83] theory of two-layer fluid flow in a inclined D . V O L A , L. B O S C A R D I N A N D J . C . L A T C H E , Laminar numerical [81] stability Phys. Fluids. (6)12, pp. 3906-3922, (1994). problem for a non-Newtonian Z. Y A N G A N D Y . C . Y O R T S O S , Asymptotic tries of large aspect ratio. solutions of miscible Phys. Fluids. (9)2, pp.286-298, (1997). 191 liquid. J . Fluid Mech. displacements in geome- Bibliography [86] Z . Y A N G , Y . C . Y O R T S O S A N D D . S A L I N , Asymptotic placements in random porous regimes in unstable miscible dis- Advances in Water Resources. (25), pp. 885-898, media. (2002). [87] S.G. Y I A N T S O S A N D B.G. HlGGINS, Linear perposed fluids. [88] C.-S. Y I H , (1967). [89] stability of plane Poiseuille flow of two su- Phys. Fluids, (31) ,pp. 3225-3238 (1988). Instability due to viscosity S.T. Z A L E S A K , Fully Multidimensional stratification. Flux-Corrected J. Fluid Mech, 27(2), pp. 337-352, Transport Algorithms for Fluids. J. Non-Newt. Fluid Mech, (100), pp 49-76, (2001). [90] E. Z E I D L E R , Applied Verlag, (1995). functional analysis [91] M. Z E Y B E K A N D Y . C . Y O R T S O S , pp. 421-442, (1992). [92] : applications Parallel flow in Hele-Shaw J. Z H A N G A N D LA. F R I G A A R D , Dispersion fluids to mathematical in a duct of large aspect ratio. cells. effects in the miscible physics. Springer- J. Fluid Mech. (244), displacement of two Accepted for publication in J. Fluid Mech, to appear (2005). [93] J.J.M. Z U I D E R W I J K , Mud Displacements in Primary .Engineers, paper number SPE 4830, (1975). 192 Cementation, Society of Petroleum
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Transient effects in oilfield cementing flows
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Transient effects in oilfield cementing flows Moyers González, Miguel Angel 2006
pdf
Page Metadata
Item Metadata
Title | Transient effects in oilfield cementing flows |
Creator |
Moyers González, Miguel Angel |
Date Issued | 2006 |
Description | In this thesis we study instabilities and mud channel removal during the primary cementing of an oil well. This process involves displacement of a sequence of non-Newtonian fluids along a narrow eccentric annulus. Previously, this has been modelled as a pseudo-steady process using a Hele-Shaw approximation. In such models, the stream function is governed by a steady elliptic problem and the fluids, (modelled via a concentration equation), simply advect along the annulus. It has been shown that for certain rheological and physical parameters, the displacement front will advances much faster on the wide side of the annulus than on the narrow side. In extreme cases the displacement front does not advance at all on the narrow side of the annulus and a static mud channel results as the finger advances up the wide side. In this thesis we consider whether the interface of a progressively advancing finger will remain stable to small perturbations. There is in fact experimental evidence that interfacial instabilities can occur in this situation. We find that the interface is in fact stable whenever there is a static mud channel on the narrow side of the annulus. Consequently, we also investigate how a mud channel might be removed by pulsation of the flow rate. Study of these two phenomena cannot be undertaken with the pseudo-steady framework. Therefore, we extend this model to flows that are fully transient. The transient model consists of a nonlinear evolution equation for the stream function. In chapter 3 we show that this transient model is in fact well-posed. In chapter 4 we study stability of multi-layer parallel flows, i.e. long fingers. If both fluids are yielded at the interface, instabilities may arise for different combinations of the 14 dimensionless parameters. These instabilities are related to a jump in tangential velocity at the interface and do not appear to have been identified before. In chapter 5 we investigate the case where a static mud channel develops on the narrow side of the annulus. Our stability theory predicts only linear stability. We therefore study the effects of a finite pulsation of the flow rates via numerical simulation. It seems that if we perturb the flow from the beginning of the displacement, the transient model fully captures the effects of the perturbation and the width of the mud channel is reduced. The pseudo-steady velocity model does not report any significant changes with respect to the results using a constant flow rate. If however, pulsation is applied after the mud channel has already formed, the removal of the mud channel will be unsuccessful. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-01-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0080123 |
URI | http://hdl.handle.net/2429/18479 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_2006-130245.pdf [ 8.46MB ]
- Metadata
- JSON: 831-1.0080123.json
- JSON-LD: 831-1.0080123-ld.json
- RDF/XML (Pretty): 831-1.0080123-rdf.xml
- RDF/JSON: 831-1.0080123-rdf.json
- Turtle: 831-1.0080123-turtle.txt
- N-Triples: 831-1.0080123-rdf-ntriples.txt
- Original Record: 831-1.0080123-source.json
- Full Text
- 831-1.0080123-fulltext.txt
- Citation
- 831-1.0080123.ris
Full Text
Cite
Citation Scheme:
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>
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-0080123/manifest