TRANSIENT EFFECTS IN OILFIELD CEMENTING 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 A b s t r a c t 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. 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 1.2 Non-Newtonian Fluids 3 1.2.1 Time independent non-Newtonian fluids \ 3 1.3 Literature Review 5 1.3.1 Industrial background 5 1.3.2 Displacement flow in different geometries 7 1.3.3 Displacements in a Hele-Shaw cell 9 1.3.4 Interfacial instabilities 13 1.4 Overview 14 Chapter 2. Model derivation 16 2.1 The pseudo-steady model and some open questions 17 2.1.1 Open problems and limitations 19 2.2 Modelling transient bulk flow cementing displacements 20 2.2.1 Non-Dimensionalisation 22 2.2.2 Reduced shear flow equations 24 2.2.3 2D averaged velocity field ; 25 2.2.4 Fluid concentration field 26 2.3 Reduced shear flow models 29 2.3.1 Stability of u and us, at (</>,£): 31 2.3.2 Stability of ua and us, at (</>, £): 33 2.4 Eccentric annular Hele-Shaw displacement transient model 35 2.4.1 Rheological assumptions 35 2.4.2 Hele-Shaw model derivation 37 2.4.3 Boundary conditions 40 2.5 Interface tracking model derivation 42 Chapter 3. Qualitative results . 47 3.1 Existence & uniqueness of ^ 48 3.1.1 Preliminary results 49 iii Table of Contents 3.1.2 Definition of dIhL2[u] and dI2^M 55 3.1.3 Steady state problem 56 3.1.4 Transient problem ••• • 57 3.1.5 Existence of a solution 59 3.2 Continuity 60 3.2.1 Preliminary results 60 3.2.2 Continuity with respect to rheological and physical properties 69 3.2.3 Steady State 70 3.2.4 Continuity for Evolution equation 76 3.3 Qualitative behaviour of solutions 79 3.3.1 End conditions and test concentrations 80 3.3.2 Steady state problem 81 3.3.3 Transient problem 84 Chapter 4. Stability of multi layer flows 90 4.1 Parallel multi-layer flows 92 4.1.1 Static mud channels 97 4.1.2 Summary of results 101 4.2 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 115 4.6 Numerical Results 119 4.6.1 Numerical Method 119 . 4.6.2 Yield stress fluids 123 4.7 Long wavelength asymptotics 127 4.7.1 Newtonian fluids 127 4.7.2 Yield stress fluids 131 4.7.3 Industrial application 138 4.8 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 145 5.1 Augmented Lagrangian Method 146 5.1.1 Principles of the method 146 5.1.2 A L G 2 148 5.1.3 Convergence of A L G 2 149 5.2 Application of A L G 2 to the displacement problem 151 5.2.1 Steady model 151 5.2.2 Transient model 154 5.3 Finite Element discretization 155 5.4 Flux Corrected Transport Algorithm 156 5.5 Numerical simulations for the annular displacement with pulsating flow rate 160 iv Table of Contents 5.5.1 Validation ...160 5.5.2 Pulsation effects 161 Chapter 6. Summary and conclusions 183 6.1 Main contributions of the thesis 183 6.2 Interfacial instabilities 184 6.3 Mud channel removal 185 Bibliography 187 v L i s t o f T a b l e s 4.1 Code validation L i s t o f F i g u r e s 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 2 1.2 Schematic behaviour of shear stress vs. shear rate for time independent non-Newtonian fluids 4 2.1 Geometry of the well 18 2.2 Example of concentration field at fixed £ and t < 4 5 3.1 The function x(|V*|,f/',7v,K;,// 0 0,n) 51 3.2 a) 7 > T T / 2 , b) (3< T T / 4 , 7 < T T / 2 , C ) R > T T / 4 , 7 < T T / 2 6 4 3.3 .0 < T T / 4 , 7 < T T / 2 65 3.4 a) Decay of ||* — ^s\\L2(n)- b) Fixed concentration field 88 3.5 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/S 89 4.1 Multilayer base flow, both fluids yielded at the interface. Rhelogical and physical para-meters 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 4 .2 Multilayer base flow, both fluids yielded 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 4 .3 Multilayer base flow, fluid 2 unyielded at the interface. Rhelogical and physical parame-ters are: 7v,i = 1, K I = 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 4 .4 Multilayer base flow, fluid l unyielded at the interface. Rhelogical and physical parame-ters 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 4 .5 Contour plot of fa>min, eccentricity vs. yield stress of fluid one, 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 .14) is fulfilled 99 4 .6 Contour plot of fa,mini eccentricity vs. yield stress of fluid one, buoyancy parameter negative. Rhelogical and physical parameters are: «i = 0 .5 , pi = 1.1, mi = 1, 7v,2 = 0 .5 , K 2 = 0 .4 , pi — 1, mi = 1 and p-= 0 . The zero level curve corresponds to the equality of condition ( 4 . 1 4 ) , i.e. faimin — 1, the rest of the level curves correspond to different values of </>i,min when condition (4 .14) is fulfilled. 99 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, p2 = 1.1, m2 = 1 and 0 — 0. The zero level curve corresponds to the equality of condition (4.14), i.e. (p^mm = 1, the rest of the level curves correspond to different values of 0i,min when condition (4.14) is fulfilled 100 4.8 Contour plot of 4>i,wini eccentricity vs. yield stress of fluid one, equal consistency. Rh-elogical and physical parameters are: K\ = 0.5, pi = 1, mi = 1, 7y ] 2 — 0.5, K2 = 0.5, P2 = 1, m2 = 1 and 0 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. 0i,min = 1, the rest of the level curves correspond to different values of <&i,min when condition (4.14) is fulfilled '. 100 4.9 Contour plot of 4>^mim 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, p2 = 1, m2 = 1 a n ^ 0 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. <fo,mjn = 1, the rest of the level curves correspond to different values of 0jiminwhen condition (4.14) is fulfilled 101 4.10 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-hand-side, 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 4.11 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, p2 = 0.75. b) pi = 1, p2 = 1.25 112 4.12 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1, P l = l, P 2 = 1 e = 0.0, & = 0.5 and 0 = 0 112 4.13 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1, e = 0.0, <fc = 0.45 and /? = 0. a) px = 1, p2 = 0.75. b) px = 1, p2 = 1.25 113 4.14 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI — 1, e = 0.0, & = 0.55 and 0 = 0. &) pi = 1, p2 = 0.75. b) pi = 1, p2 = 1-25 113 4.15 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, p2 = 0.75. b) pi = 1, p2 = 1.25 114 4.16 Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1.5, K2 = 1, e = 0.0 and /? = 0, mi = 1.5, m2 = 2. a) = 0.5. b) <pi = 0.6. For b < 0, pi = 1 and p2 = 0.75, for b = 0, pi = 1 and p2 = 1 and for b > 0, pi = 1 and p2 = 1.25. .. 117 4.17 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, m2 = 2. For 6 < 0, pi = 1 and p2 = 0.75, for b = 0, pi = 1 and p2 = 1 and for 6 > 0, pi = 1 and p2 = 1.25. . . . . . . . 118 4.18 Maximum imaginary part of s vs. a. Inclined annulus. Rheological and physical para-meters are: KI = 1.5, K 2 = 1, mi = 1, m2 = 1, b = —1 and 0 = 0. a) Exact solution, b) Numerical solution 121 4.19 Maximum imaginary part of s vs. a. Rheological and physical parameters are: «i = 1.5, « 2 = 1, e = 0.0 a n d 0 = 0, mi = 1, m2 = 1, ry^ = 0.9 and 7y, 2 = 0.7. a) (pi = 0.5. b) (pi = 0.6. For 6 < 0, pi = 1 and p2 = 0.75, for b = 0, pi = 1 and p2 = 1 and for b > 0, pi = 1 and p2 = 1.25 123 viii 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, m2 = 1, TY,I = 0.7 and 7V,2 = 0.9. b) K\ — 1.5, K 2 = 1, e = 0.0, fa = 0.5, 3 = IT/8, mi = 1, m2 = 1, 7v,i = 0.7 and 7v ,2 = 0.9., mi = 2, m,2 = 1-5. For b < 0, pi = 1 and p2 = 0.75, for 6 = 0, p\ = 1 and P2 = 1 and for 6 > 0, pi = 1 and p2 = 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 2 = 1> e = 0.0 and /3 = 0, mi = 1, m2 = 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 p2 = 0.75, for b = 0, pi = 1 and p2 = 1 and for b > 0, pi = 1 and p2 = 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, m2 = 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 / 2 , 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 / 2 , dashed line shows the imaginary part of fe and fe and dashed thick line shows the magnitude of fi and fe 126 4.23 Contour plot of s2. Rheological and physical parameters are: K\ = 1, pi = 1, e = 0.0, fa = 0.5 and P = 0 131 4.24 Contour plot of s2, «i vs. K 2 . Rheological and physical parameters are: pi = 1, mi = 1, m2 = 1, 7v,i — 1, TV,2 = 1» P2 = 1, e = 0.2, fa = 0.5 and P = 0 136 4.25 Contour plot of s2, K\ VS. K 2- Rheological and physical parameters are: mi = 1, m2 = 1> 7v,i = 1, 7v ,2 = 1, e = 0.2, fa = 0.5 and P = 0. a)pi = 1.1 and p2 = 1. b)pi = 1 and p2'= 1.1 137 4.26 Contour plot of s2, n\ VS. K 2 . Rheological and physical parameters are: pi = 1, TV,I = 1, 7v ,2 = 1, P2 = 1-1, e = 0.2, (pi = 0.5 and P = 0. a) mi = 1.2 and m2 = 1. b) mi = 1 and m2 = 1.2, 137 4.27 Contour plot of s2, K,\ VS. K 2 , buoyancy parameter equal zero. Rheological and physical parameters are: KI = 1, pi = 1, /t2 = 1, p2 = 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 2 : (a) « 2 — 0.5; (b) re2 = 1; (c) K 2 = 1.5. Fixed parameters are: TV,I = 1-0, K.\ — 1.0, mi = 1, m2 = 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, m2 = 2, b = -0.5, fa = 0.5 and P = 0. a) K 2 = 0.5. b) K 2 — 1- c) K 2 = 1-5 141 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 2 = 0.4, mi = 1, m2 = 1.2, pi = 1, p2 = 0.9, TV,I = 1, TV.i = 0.9, TV,I = 0.7 e = 0.0 and P = 0 161 5.2 Initial condition for the displacements 162 ix List of Figures 5.3 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2-K/UI, U = 10 and Sp = 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, m2 = 1.2, p\ = 1, p2 = 0.9, 7V,i = 1, T V , I = 0.9, 7y,i = 0.7 e = 0.3, 0 = 0 165 5.4 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Steady state velocity model. Period T = 2ix/u), u> = 10 and Sp = 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, K2 = 0.4, mi = 1, m2 = 1.2, pi = 1, p2 = 0.9, T V , 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, 5P = 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, K2 = 0.4, mi = 1, m2 = 1.2, pi = 1, p2 = 0.9, ry,i = 1, T V , I = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 167 5.6 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2-K/LOC, to = 10, Sp = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 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, 5P = 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, K2 = 0.4, mi = 1, m2 = 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 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2w/uje, w = 1, Sp = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, ry,i = 1, T V , I = 0.9, T V , I = 0.7 e = 0.3, 0 = 0 170 5.9 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 27r/we, UJ = 0.1, 5p = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. Physical and rheological parameters: KI = 0.5, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 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 5P = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 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 5P = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, 7v,i = 1, T V , I = 0.9, 7Y, i = 0.7 e = 0.5, 0 = 0 173 x List of Figures 5.12 Displacement flow in eccentric annulus, interface propagation. Unsteady displacement. Transient velocity model. Period T = 2rr/W, w = 10, 5P = 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: « i — 0.5, K 2 = 0.4, mi = 1, m2 — 1.2, p i = 1, p 2 = 0.9, 7V,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 8P = 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 2 = 0.4, mi = 1, m2 = 1.2, p i = 1, p 2 = 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 Sp — 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 2 = 0.4, mi = 1, m2 = 1.2, p i = 1, p 2 = 0.9, 7 v , 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 5p = 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: « i = 0.5, K 2 = 0.4, mi = 1, m2 = 1.2, p i = 1, P2 = 0.9, T V , 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, 6P = 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, K2 = 0.4, mi = 1, m2 = 1.2, p i = 1, p 2 = 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, 5p = 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, « 2 — 0.4, mi = 1, m2 = 1.2, p i = 1, p 2 = 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 5p = 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, /t2 = 0.4, mi = 1, m2 = 1.2, p i = 1, p 2 = 0.9, T V , I = 1, T V , 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, 5p = 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, K 2 = 0.4, mi = 1, m2 = 1.2, p i = 1, p 2 = 0.9, T V , I = 1, T v i = 0.9, T V I = 0.7 e = 0.8, 8 = 0 ' 182 xi A c k n o w l e d g e m e n t s 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 collab-orative 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 I n t r o d u c t i o n 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, spacer fluids and 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='ni crxz=ayz = 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. Any 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 non-Newtonian 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: (1.2) V-Voo Vo -Vco (1 + [Ki)m) which is called the Cross model, and V-Voo 1 (1.3) r?o -Voo (1 + (tf 7) 2)m / 2 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 (1.4) 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, ri = p + j (1.6) 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 = v/7y"+ s/K^ (1.7) which is called the Casson model, and is given in terms of the shear stress, and finally r, = KT~l + ^ (1-8) 7 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 Literature Review 1.3.1 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 under-stand 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 divi-sion 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 ap-proximated 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 to find the 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. The fluids considered satisfy the Herschel-Bulkley model. Because the annulus geometry is long and thin they derive a 2-D model of the bulk fluid motions 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 equation field will determine the position of the interface be-tween the fluids. 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 long finger develops 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 be found analytically for concentric annuli and annuli with small eccentricities via a perturbation method. Their ap-proach 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, for when 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 displace-ment 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 two fixed parallel plates which are sufficiently close together (Hele-Shaw cell) is u = -^-^-zld- z) and v = -^-^-z(d - z) (1.9) 2pdx y 1 2pdy y ' K ' 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 law for a 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 of flows. 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 thinner fluid penetrates the thicker fluid in 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 of fingers. The displacement is ineffective and the interface is said to be unstable. This phenomenon has been termed "viscous fingering". 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 the fluids is, 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 this field of work involves both immsicible and miscible displacements, we concentrate in the latter ones which are the type of flows considered in the thesis. In 1986 Tan & Homsy, [74] studied the linear stability of a miscible displacement in a porous medium for rectilinear flow. They assume that the fluids have 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 displacing fluid is 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 a finger becomes large enough the tip of the finger becomes 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 = -K1-^)^ | vp |>G' (i-n) « = 0 |Vp| < G (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 can be found i n the l i terature, see for example [4, 16, 34, 44, 53]. Ac tua l l y , the mode l derived i n [10], wh ich is the model we follow and extend i n this thesis, uses a nonlinear Da rcy ' s law w i t h a l i m i t i n g pressure gradient as wel l . H a v i n g stated this, we review some of the work that has been done concerning non-Newtonian fluid displacements i n a Hele-Shaw cell or a porous med ium. W i l s o n , [84], considered the Saffman-Taylor problem for a power law f luid displaced by air. B y means of linear s tabi l i ty analysis he found that the growth rate of the ins tab i l i ty is 0 ( n - 1 / 2 ) , where n is the power law index. La t e r Aza iez & Singh, [4] considered miscible displacements of shear th inn ing fluids by a Newton ian fluid or vice versa. T h e i r approach is the same as the one reported by T a n & Homsy, [74] for Newton ian fluids, the quasi-steady-state approach. T h e y found that for the n N - N case, (meaning non-Newtonian-Newtonian displacement) , the displacement can be unstable even for favorable mobi l i ty rat io, this was expected due to the shear th inn ing behaviour of the f luid. For the N - n N case they found s imi lar results as for the N - N case. Studies i n wh ich bo th fluids are non-Newtonian are presented below. Pasca l s tudied the classical planar interfacial instabil i t ies i n porous media for immisc ib le visco-plast ic fluids, [53, 54, 55]. He used the Muska t approach, i n which one assumes a long finger that extends ahead of the moving interface, i f the veloci ty of the fluid w i t h i n the finger is faster than the veloci ty of the fluid outside, the interface is unstable. B y doing this he showed that i n contrast w i t h what 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 an unfavorable mobi l i ty ratio, i.e. less viscous f luid displacing more viscous fluid, a posit ive difference of y ie ld stresses w i l l a id the displacement. B y positive we mean that the displacing fluid should have a higher y ie ld stress than the displaced one for this to happen. H e also considered the effects of inc l ina t ion of the media , and showed that gravi ty forces tend to stabil ize the interface as long as the heavier f luid is under the l ighter one. In 1999 Coussot , [16], considered the Saffman-Taylor problem invo lv ing two y ie ld stress fluids. H e found the same characteristics of instabil i t ies as for Newton ian fluids except that the wavelength of m a x i m u m growth can be smal l even at vanishing velocities, i.e. ins tabi l i ty may occur even i n the absence of gravi ty effects, as soon as the y ie ld stress of the displaced f luid 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 half-plane. Yiantsios and Higgins [87] consider a plane Poiseuille flow, in which 2 fluids flow together 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 shear flow of two superposed miscible fluids with 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] found similar results using a different profile for the viscosity variation. For the case of multi-layer flows in porous media, Raghavan & Marsden, [62] showed, in the presence of density difference, if the heavier fluid was in top of the lighter one, the flow was 13 Chapter 1. Introduction 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). 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 immiscible fluids in 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 the fluids. 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 parallel flow of a gas on top of a viscous fluid in 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 geome-try 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 the fluids is 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 parallel flows in 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 displacement flows presented 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 displacement flow and 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 shear flow model from the Navier-Stokes equations and discuss a number of the modelling assumptions. In §2.3 we approximate the reduced shear flow model from §2.2 by a gap-averaged model. We show that this approximation is valid when the concentration field changes slowly. In §2.4 we finally derive the transient Hele-Shaw model from the gap-averaged reduced shear flow model 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 concentration field. This variant is more convenient for studying interfacial stability. Since the concentration field is 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 in 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 in [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-[Sa + f] = 0, (2.1) X(|V*8|) T V W s <=> \Ss\>jj, (2.2) | V * 5 | = 0 <=* \SS\<^, (2.3) where \&s is the stream function. The unwrapped narrow annular space is (0,0 € (0,1) x (0, Z), the annular gap half-width is H{<j>) = 1 + ecos7T0; e G [0,1) is the eccentricity, defined as A r / ( r 0 — ri), see Figure 2.1 . The vectorfield / and the rheological properties of the fluids may 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 / S | , arising from the viscous shear-thinning behaviour of the fluids, and Ty is the fluid yield stress. The concentration is simply advected: l[Hc} + -^[Hvsc] + ^-z[HWsc} = 0, (2.4) - W = H w s , -j; = -Hva. (2.5) 17 Chapter 2. Model derivation Figure 2.1: Geomet ry of the well T h e der ivat ion of (2.1), at fixed t, follows straightforwardly from the reduced shear flow model : 0 = - ^ + — T a m d m + g 6 , (2.6) (2.7) dp d _ dp d where r 5 = (Ts^y,TS£y) depends on the veloci ty gradients (vSty,wS:y) at leading order, and (<?<£> 5£) represents the gravi ta t ional acceleration. These equations, plus the incompress ibi l i ty condi t ion , are essentially averaged w i t h respect to y, that is, the radial coordinate scaled relative to the distance from the centreline of the annulus, and cross-differentiated,-to give (2.1). T h e pseudo-steady model above has been studied i n some deta i l . In a uni form inc l ined eccentric annulus, i t appears that the displacement is characterised by one of 3 different behaviours. a) If the interface between fluids advances steadily along the well at the same speed a l l a round the annulus, this leads to a good displacement. These steady states have been observed i n the simulations of [10]. For smal l eccentricities, ana ly t ica l expressions can be found that describe their shape, see [57]. T h e i r s tabi l i ty has been characterised i n [58]. b) If a steady state displacement is not found, then invar iably the interface remains unsteady and elongates into a progressively long 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 displacement flows are 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 r0, respectively. fV is just the radius of the inner cylinder and f0 is the radius of the outer cylinder, see Figure 2.1. We define the mean radius as f* = \(f0 + fj) and the mean half-gap width d = \{r0 — fj). 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 „ A 19 . . . . 1 d „ d „ rgg dp , n „ . \dt J rdr r 06 d£ n r or Jdv . \ 1 a , . 2 . , 1 9 . 0 . ldp _ .... P ^ + u - V v ) = - - { r \ e f ) + - - r 6 6 + - , M - - f - ^ p ^ (2.9) Jdw „ „ \ l f l . . . , 1 9 . dp _ • 0 = ^ du, rdry ' rd9 d£ where the gravitational vector is given by g? = gsm(3{€)cos6, gg = g sin f3(£) sin 9, g^ = gcosf3(£). | + + + = V . ( 6 ( c « ) V c ) (2.12) where c € [0,1]. 21 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 = —, <P = -irr* 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[ffc,y + kk(Y)nk], p.* = Z-, p* = ^f. k j* 0 The dimensionless rheological parameters are defined by: and the fluid densities are scaled with 23 Chapter 2. Model derivation p*=max[pk]. 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*af 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 / (2-13) dy dv 1 dp d , . pm = -7ad4> + d-yT^ + sT*' ( 2 1 4 ) dw dp d ge . pm = -i + d-yT^ + s%' (2A5) dy ra d<p <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: u = v = w = 0, at y = H, (2.17) u = T(f>y = T^y =0, at y = 0, (2.18) 24 Chapter 2. Model derivation The leading order rate of strain is characterized by: 7 (2.19) In regions where the fluids are yielded we have: T4m V(i) + 0(5*M riy ~ V(j) + 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 vdy wdy. 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" 10m2s"\ thus Pe » 1 and typically Pe ~ 1010 -1012. Denoting the scaled 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*v/t*av where t*av = jff by e, i.e. ._ p*5*2irr*aU* A* Formally e —> 0 as 5* —> 0. However, although t* = 0(S*2), the ratio t*/t*av is often 0(1). This 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 ~ 108 — 1010 and the aspect ratio in a typical well can not justify ePec <S 1. 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 PeG1^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 one fluid by a lighter and less viscous fluid, between parallel plates. If frontal instabilities don't occur and a finger develops, there have been many studies of interfacial instabilities in shear flows, see e.g. [5, 42, 67] for an overview. In general, such flows are not robustly stable to all wavelengths and rheology combinations. The situation is not apparently helped by having a thin diffuse interface layer over which fluid properties change. Recently in, [23] the authors studied miscible multilayer flows where 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] and find a transient dis-placement 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 approx-imation 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 = G + ^ - r , (2.25) dy where u = (v, w) and G = (Gh GS) EE + g$, - | + g^j . (2.26) The system (2.25) is supplemented with the conditions (2.13) & (2.16) and (psin/3sin7T0 pcosjA (9*>9t)={ sf (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 = Gs + — T S . (2.28) dy The subscript s is used to indicate that we are seeking a velocity field us = (vs,ws) that satisfies (2.28), where potentially GS ^ G, (e.g. because the fluid concentration will be advected differently by us than by u). As with (2.25), the pressure and concentration are independent 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): . pua,t = Ga + ^ T A , (2.29) where again the subscript a denotes that we seek a velocity field ua = (va,wa), and that potentially Ga ^ G; (2.29) is supplemented with the incompressibility condition: dua dva dwa The overbar indicates a gap-averaged quantity, e.g. ua = (va,wa): V a = H J0 V a y ' W a = H j W 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 G3 at a point ((f), £), we are able to find a unique solution (vs,ws), (and could then integrate the incompressibility condition to find ua). All the time dependency in (vs,ws) must come from the vectorfield Gs, which is determined from (2.1), with appropriate boundary conditions, and the advection of the concentration field. When Gs and cs vary only slowly with time, (vs, ws) are (pseudo-)steady on the 0(1) timescale. 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 (va,wa), but it is slightly unclear if, for specified Ga and ca, the evolution of (va,wa) is properly determined. Rather, given time derivatives of the averaged velocities (vatt,wa,t), w e c a n determine unique (va,wa), (i.e. this is the same problem as (2.28) 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_1. It will however be the gap-averaged system (2.29) that forms the basis 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 u8, at (d>, £): For n = 1,'problems (2.25) and (2.28) are of a form considered by [20, 30]. Multiplying (2.28) by u.s and integrating with respect to y, we have: as(us, us) + js(us) = L(Gs-us) (2.31) where s(u,v) = KS / Uy-Vy dy, Jo (2.32) From this we can deduce that: f" js(v) = TY,S / \vy\ dy, (2.33) Jo L{f) = [ fdy. (2.34) Jo \us\\m<—\\Gs\\L,, ' (2.35) KSTT from the Cauchy-Schwarz and Poincare inequalities, where | |M s | | i / i is the semi-norm: 1/2 Jo Mm = ( / N 2 dy which is equivalent to the H1-norm. Using the methods in [20, 30] we may assert that there exists a unique solution us. Time-dependency enters us only through Gs and cs. We assume that us is differentiable with respect to each of the physical fluid properties, which depend on cs, and with respect to the components of Gs. Therefore, we may suppose that uStt = eb(et), (2.36) for some vector b. 31 Chapter 2. Model derivation Now we consider stability of (2.25). Multiplying (2.25) by (us — u) and integrating with respect to y we have: p ut- {us —u)dy > L(G.(us — u)) — a(u, us — u) Jo -j(Us)+j(u), (2.37) where a and j are defined analogously to as and jg, but with non-subscripted rheological constants. Combining (2.36) with (2.28), we have: puSft = Gs + Ts,y + epb. (2.38) Multiplying by (u — us) and integrating with respect to y gives the following inequality: fH p / uSit • {u - u3) dy > L(GS • (u - us)) - as(us, u - us) (2.39) Jo -js(u) + js(us) + epL(b • (u - us)). Summing the two inequalities we have o"57 / \ u s ~ u \ 2 & y < -a{us - u,us - u) + L((G - Gs) • {u - us)) z at J0 +epL(b • (us - it)) + ^1 - -^j as(us, u - us) + ( l - ^ ) [ ? - ( « ) ( 2 - 4 0 ) Using the Cauchy-Schwarz inequality we have: | ^ | | U s - « | | 2 a < - K \ \ u , - u \ \ 2 H i (2.41) + (\\G-Gs\\L2 + ep\\b\\L2)\\us-u\\L2 + (\KS - K\ \\US\\HI + \TY,S -TY\H1/2^ \\US - u\\Hi. We may split the first term on the right hand side in two. From the Poincare inequality: \\us-u\\Hi>^\\us-u\\L2. Therefore either we have that: 2H (\KS-K\ \TY,s-rY\\ / 0 4 0 x I K - < — \-^r\\G'\\» + gi/2 ) . ( 2 - 4 2 ) 32 Chapter 2. Model derivation or we have the decay equation: d K7T2 p— \\us -u\\L2 < - ^ \ \ u s - u \ \ L 2 + (\\G-Gs\\L2+ep\\b\\L2), from which, for t € [0, T] \G-GsWv |« s-u||L2(t) < / Jo + ep||b||L2 e-Xi(t-s) d g +||t t s-n||L 2(0)e-A l t, (2.43) with A r = m i n ^ ^ J . • f ^2 ) lin te[o,T] In either case, we note that locally 1\us — M|| jji is bounded by the difference in the data between steady and transient problems. 2.3.2 Stability of ua and ua, at (tp, £): We proceed analogously to the above, with the variational form of (2.29) & (2.28). The latter we modify, to include the term us<t, i.e. pHus,t • (ua - its) > L(GS • (ua - us)) - as(us, ua - us) (2.44) -js(ua) +js(us) + epHb • (ua - us). pHua,t • (us - ua) > L(Ga • (us - un)) - aa(ua, us - ua) -ja(us)+ja(ua). (2.45) Summing these two inequalities, we have * 2 ^ ^ l l M s - « o | | i i < - a a ( u s - u a , u s - u a ) + epHb-(us-ua) +L((GS - Ga) • (us - ua)) + ^ l - ^ j a s ( u s , u a - U s ) + ( l - ^ j [ j s ( u a ) - j s ( u s ) } . (2.46) 33 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\\GS\\L2 . .2H1'2 \\Us-Ua\\Hi < \Ka-Ka\ h \TY,s ~ TY,a\ 7TKsKa Ka + \\Ga-Ga\\v> — , (2.47) 7TK a or we have the decay equation: p d ,. KaTT2 . r from which, for t G [0, T] | | « 8 - « a | | L i ( t ) < f eH\b\e-x^-^ As Jo ' + I K - « a l b ( 0 ) e - A 2 t , (2.48) with , . f HA7T2 \ A2 = mm i 0 }. t€[0,T] \2H2pj Again we see that ||u s — w||x,i remains uniformly bounded. Note that we cannot expect results on \ \us — ua\\L2, since the averaged problem leads only to an equation for evolution of the mean velocity. Summary and further comments Combining (2.43) & (2.48), we see that ||tt0 — u||Li(t) 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 | |« a — ,u||ii(t) = 0(e). Since we 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-1, on which timescale c also evolves. Decay of initial errors occurs on an 0(1) 34 Chapter 2. Model derivation timescale, after which the errors are simply due to the differences in G and the material prop-erties. These are governed by the evolution of the gap-averaged velocity fields, which takes place on the slow timescale. We should therefore expect \\ua — u\\Li(t) = O(e) on a timescale [0,T], where T~l ~ o(e), but T - 1 ~ 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 \\ua — = 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 2.4 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 ~ wy Thus, \\u — v\\Hi terms and the Poincare inequality on W1,n+1(0,H) for the rest. 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, (2-49)' 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: 7 i j . T > T V , 7Vi — 1 (2.50) 7 = 0, 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 . lij, T>Ty, (2.51) 7 = 0, r < 7 v , 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 high-shear viscosity we are later forced to work in subspaces of L1+n and Wl,1+n, rather than I? and H1. Although most of what we prove in the Hilbert space setting can be extended to L1+n and Wl'1+n, it is less convenient and involves additional analysis. From the numerical perspective this is anyway unnecessary, since numerically we end up working in finite dimensional subspaces, which lie in both {L1+n, W^1+n} and in {L2, H1}. 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 static fluid regions in the annulus is commonplace in cementing. Fluid concentration effects Finally, if we have intermediate fluid concentrations occurring in the annulus, it will be neces-sary 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 prop-erties. Whilst this is certainly not realistic for all fluids combinations, 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) = T Y I 2 C + 7-^1(1 - c), K(C) = K2C + KI(1 - c), (2.52) n(c) = n2c + ni(l — c) p(c) = P2C + pi(l -c). (2.53) 2.4.2 Hele-Shaw model derivation We first integrate (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) = [Ga-pU, a,t] I y v(y) ^ dy, where u = (v, w) and Ga dp psmBsinircf) dp pcosB "d~4> + St* '~"di St*~ (2.55) (2.56) Therefore ua is instantaneously parallel to Ga - pua,t- Writing s(y) = \ua\{y) and A = \Ga — pua,t\, at each fixed (<fi,£,t), the speed s is related to the modified pressure gradient A via the one-dimensional shear flow problem: A dy ds dy 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 r-H j r-AH = JQ yi{y) dy = ^2j T^(T) d - r (2.57) where from (2.51), for r > ry, and 7 ( T ) is obtained by inverting this monotone function. We may observe straightforwardly that 7(r) is a strictly monotone C°° function of r > ry. If n < 1, then as T —> ry, 7(r) —> 0 and the Herschel-Bulkley term dominates: i 1/n T — Ty * 1 7 As r —> oo, the high-shear viscous term dominates: T — TY 7 Mo In between these limits numerical inversion is needed to define J(T). If n = 1, then T — Ty 7 : Poo & V r > r y . 38 Chapter 2. Model derivation The expression (2.57) defines the closure relationship between the gap-averaged flow rate | 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. | W | = I W K x j t f . T Y . K . n . / Z o o ) , x = x(\W\;H,Ty,K,n,p,00). Returning to (2.55), we average across the half-gap, y £ [0, if], to give: *. = [Ga-puaA±[ j" -faiyty, (2.59) from which we see that ua is also instantaneously parallel to the vector [G — pua,t]' <9# / <92# dp d2V dp dz'd<t>) \pd£dt e<i>+94" pd4>dt d$+9i | W | A For | V * | > 0, replacing A with x + ry/H, this implies that <92* dp psmfSsmiTcj) _ (x(|V$|) + ry/H\ <9# (2.60) d£dt d<p St* \ | V * | J di _pd2V dp pcos/3 _ [x(\V*\) + TY/H\dV d(f>dt d^ St* V J d(f> 39 Chapter 2. Model derivation Cross-differentiating to eliminate the pressure, we finally arrive at: Yx(|v*|) + 7 Y / f f V • [pV*t] - - v • /here LV iv*| v * + f (2.61) / . scos/3 . , sin/5"sin7Tc/>\ ,„ „ „ x f=(p( -c )^^)-V^j- (2-62) If | V*| = 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 |S|>,y/2f, (2.64) |V*|=0 \S\<TY/H. (2.65) We note also that: d24> dp pcosB d2^ dp p sin B sin nq depot <9£ St" ' ^<9£(% <9c/> S i * (2.66) Equation (2.63) is the classical formulation of the evolution problem for (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. (2.67) 40 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$> • S5 = O^^(0 ,O,i)=O. (2.69) a* 55 = 0^— (<f,,Z,t)=0. (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 * 0 ut( < A) at z = Z, respectively. We might then impose *(>,(),*) = * in(0,t), (2.72)' tt^.Z.i) = *„*(>,*), (2.73) 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 two fluids of different constant physical properties. Below we derive this model. Derivation and jump conditions at the interface From §2.4.2 we have, V-[pW t] = - V •[£ + /] (2.74) Let us assume for the moment that * e C1(il). Now, let us multiply (2.74) by a test function v such that, v : R2 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 = / ( P V*i + S + f)- Vvdfl = - / V • (pVtf* + S + f)vdfl (2.75) Jo. JQ We derived this equality assuming that * £ C1(fi), but note that (2.75) has meaning even if 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 only fluid 1 on it and 02 that part on the right of the curve with fluid 2 on it. Suppose that ^ is a solution of (2.74), as we will prove in chapter 3, and assuming that its first derivatives are uniformly continuous in fl\ and 42 Chapter 2. Model derivation fl2 we can choose a test function v with compact support in fii, thus (2.75) becomes, 0 = / (pWt + Si + fi) • Vvdfl = - / V • (pVtf* + Si + fi)vdfl, (2.76) JUi Jili where the subscript 1 means that we are considering the properties of fluid 1 only. Because fli is full of fluid 1 then V • fi = 0, thus PiA$t =-V • Si in fii. (2.77) Likewise, p2At>t = - V • S2 mfl2 (2.78) 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*t + S +/) • Vvdfl = - f (piVtft + Si + /1 ) • Vudfii Jn Jni + [ (p2V*t + S2 + /2) • \7vdfl2. (2.79) Jih Now, since v has compact support within fl, we have 0 = / (piW t + Si + fi) • Vvdfl! = - / (piA*t + V • (Si + fi))vdflx Jn Jni + f {(piWt + Si + fJ-n^vdl. Jc = [ {(pi^t + Si + fi)-ni}vdl (2.80) Jc because of (2.77) and ni denotes the outward normal to C from fli to f22. Similarly we have 0 = / (p2V*t + S2 + f2) • Vudft2 = - / (p2A*t + V • (S2 + /2))i;dr2x 7n Jfi2 + / [(p2Wt + S 2 + /2) • n2}vdl Jc = / [(p2V*t + S 2 + /2) • nzjvdi (2.81) Jc Adding (2.80) and (2.81) we have 43 Chapter 2. Model derivation (pfcVtft + Sk + fk) • n\\ = 0, (2.82) along C. Thus, our interface model can be written as follow, p f eA* t + V-S f c = 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 d*dF WdF 7 a r r - ^ a 0 + 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, pfc os/3 pfesin/3sin71-0 , c *|2 = 0, (2.88) VF| 2 = 0. (2.89) The first of these is simply continuity of 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 inter-face 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/Pe1/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 t), we have that c(<f>, £, t) = 1 for 4> G [0, <S>i{(,,t)) and V therefore 45 Chapter 2. Model derivation Jo. and ail ^ w = y 0 -^r^m-^ (2-91) at <{> = 0 M , (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 n ± , ,„ - ^ + ^ 7 + 7 ^ : ^ F = 0 at = 2.96 e OT oi o<p oi Therefore a natural choice for F(<j>, £, i) defined in (2.86) will be, F{fai,t) = (j)-fa{i,t). 46 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);fr1(fi))> * t €Loo([0,co);ir1("))-ii) In §3.2 we consider the flow of a fluid with 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 density field P2, 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, then for t G [0,T] | | p i - P2||i«»T + | | V ( * i - * 2)||L*(0) > ^ 2 | | V ( * i - tt2)||£9 (3.1) 2. If 7 y i —> 7 y 2 , we have liny - r 2 y | | L 2T+ ||V(*i - *2)||L*(0) > C 2||V(*i - * 2 ) H L 2 (3-2) 3. Finally if K\ —»/c2 we get I K - « 2 l U a r + | | v ( * i - * 2 ) H L * ( O ) > c2||v(*i - * 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 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. 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(* - * S ) | | L 2 < ||V(* - *s)\\i?Q)e-Kt where 1C K = for C<1. Pmin 3.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 The flow rate 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 0 , , dp pcos/3 dp psin/Jsmmft S + pV* t = ( _ _ _ ^ _ , _ _ ) (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 con-fined to the interior, at times we shall also assume the following. A4 As £ —> 0 and £ —> Z, the physical properties of the fluid, TY, p^, K, p and n all approach constant values. 3.1.1 Preliminary results We commence with a number of preliminary results. Proper t ies of %(|V*|) Propos i t ion 1 For |V*| > 0 we have that x(|V*|) is C°°, strictly positive and strictly monotone; x(|V*|) -> 0 as |V*| -» 0. Proof This follows from the properties of A(|V*|). See the discussion in §2.4.2. 49 Chapter 3. Qualitative results Proposition 2 The function x(|V*|) is bounded below by XN(\^^\), XB(|V\I>|) andxHB{\^^\), defined implicitly for |V\IJ| > 0 as follows: | V * | = Z?*L, ( 3 - 5 ) ' 3Moo g 3 X / j ( X B + 1 . 5 7 y / f f ) 3/ioo ( X B + r y / F ) 2 |V*| = J ; . A B ^ ' / _ T / , / . (3-6) V * = 7 — " T O N - 7 _ f , m 2 U g j 3 + ~ T T ^ ' m = l / n (3 .7 ) P r o o f These results follow from (2 .58) and the definition of 7(1"). We prove only ( 3 . 6 ) , the others following in analogous fashion. We bound T ( 7 ) below for r > ry, by neglecting the term KJN, i.e. T — Ty T{J) > Pool + ry < . Poo Inserting into ( 2 .58 ) : |V¥| < I [ X H + T Y T(T - ry) dr = ^ i ^ . Therefore, V|V\P| > 0 we have XKXB + 1 . 5 7 y / g ) X 2 ( X + l - 5 7 y / g ) ( X B + T y / F ) 2 - ( 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 7 N . The bound 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*| n , whereas evidently XN ~ |V*| and X B ~ I V\I/|. 5 0 Chapter 3. Qualitative results Figure 3.1: The function x( |V*| , i f ,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 n } + TV, implies that j(r) > 7m ("?"), where 7m (T) = min • T — Ty 2/ioo r — r y 2K l / n " Therefore, defining XM implicitly via | V * | = •TY/H]2L XMH+TY Tjmir) dr, [XM + r / " j J ? y leads to x(|V*|) < XM( |V* | ) . Further, we can see that at large r, T — ry 7m (r) 2/Joo that g 3 X 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 Behav iour 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 d n = / *[pV* t + S ] - i / d s - / V * - S - * V - f dfi, Jn Jen Jsi (3.8) where u denotes the outward normal to Vt. Using (3.4), we have / *[pV*t + S] v ds = / * ( - ^ + ^ , ^ - ^ ) . i / d s Jdii Jan V o<P ) < \\Vp-(9<t>,9s)\\L~(dn) I 1*1 ds JdQ < C o | | V p - ( 5 ^ f l € ) | | i o o ( 8 n ) | | * | | f f j -(3.9) The first bound follows from assumptions Al & A2, using the Holder inequality. The last line follows from, fL <9* / ffidn > / |* i n |d0+ / |* o«t |d0 / ffidfi > / |*(U,i) |d£ + / |*(0,£,i)|d£ =» 2 / |V*|dft > / |*|ds JO, Jdfl From proposition 2 we have that X(|V*|) > ^ | V * | . (3.10) 52 Chapter 3. Qualitative results Combining all this: , A ^ J | V * | 2 d f t < Co||Vp-^,5€)||icc (8n)||*||i/i + | |V-f| | L 2 | |* | | L a -- 3 ig F{^}ll*llia- is f{^}ll^J. where 1*1 HI = 7 JQ IV*I 2 dfl 1/2 Since * = 0 along c/> = 0, the seminorm ||*||#i is equivalent to | | * | | f f i , and evidently 11*| 11,2 < | | * | | / f i . Therefore, we can find constants C\ > 0 and C 2 > 0, for which d di ; j J | V * | 2 d f t < Ci | |* | | H i - C 2 | | * | | ^ i . (3.12) Integrating (3.12) with respect to i, we can find C3 > 0 for which C3\\n%i(t) < inf{f} | |* | |^( t ) JQ Z t=o Jo We see that the integrand becomes negative if l * I M * ) > Ci c2 and consequently ||*||#i(i) is bounded for all t > 0. • I(u) and dl[u] We denote by Vj the subspace of Hl(Q) containing functions that satisfy boundary conditions (2.67), (2.68), (2.72) and (2.73). The space Vj is non-empty, since for example ** G Vj, where ,T,* ,T, IJ. *\ P0Ut(<P) - P , ,Tt. ,x ^ P~ Pin{4>) / n i m * = #in(</>,i) — —- + yout(4>,t)- j— TTT, (3.13) PoutW ~ Pin(<P) PoutW) - Pin{<P) and Pin(<f>) & Pout(4>)y a r e the density at the inflow and outflow, £ = 0, Z, respectively. Note that the boundary streamfunctions, & *out, must both satisfy * i n(l,i) = * o u i(l,i) = Q(i), 53 Chapter 3. Qualitative results *i„(0,t) = *out(0 )t) = 0. We also denote by V/o the subspace of i / 1 (fi) containing functions that are zero at dfl and equipped with the L2(fl) inner product. Note that Vip is a Hilbert space and in fact Vj^ = HQ(Q), but for notation we work with Vifi. Obviously Vj is an affine space, i.e. Vi = V + VIfi, for any £ Vf. For <J>* G V/ and u € V^ o we consider the following functional, I(u): I[u] := { hlu] + h[u] u € Vi,o, (3-14) where with Ik[u] = f Lk(Vu) dfi, k = 1, 2, (3.15) i l ( V u ) = 21 ^ ^ d s + V^'+^-f , (3.16) L2(V«) .= -^|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 H0\n) : I[w] > I[u] + {v, w - u) V w G V/,0}, (3.18) where (•, •) is the inner product in L2(tt). Proposition 3 The functional I(u) is strictly convex, proper and lower semi-continuous; dl[u] is monotone. Proof 54 Chapter 3. Qualitative results For convexity we may follow the same steps as in [57], (proposition 2). Note that L(Vu) = Li(Vu) + L2(X7u) is convex. Therefore, (see [17]: theorem 3.1, p. 66 and theorem 3.4, p. 74), L(Vw) is weakly lower semi-continuous. The lower semi-continuity 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 dI2,L2[u] We now consider characterization 1 of dl[v]ij2. We consider u,w € H$(Cl). The functional h[u] is Gateaux-differentiable. Consequently, 91^2[u] = I[ L2[u], and we find that / i M > h[u] - Jn(w - u ) V • ( > j ^ J ( v * * + V«) + f ) dil (3.19) for u G V]:o, \/w G Vjfl-If u G Vifi, then dl2,L2M *s characterized by: (ry (V** + VM) ' \~H |V#* + Vu| (see [1]), which is set-valued for |V** + Vu| = 0. Thus (3.20) We now combine these 2 terms. Propos i t ion 4 Let I[u] = Ii[u] + I2[u], with h[u] Gateaux-differentiable and dli[u] single valued. Then if v\ G dKi[u] v = v\ + v2 G d(h[u] + I2[u\) & v2 G dl2[u). Proof JThe sub script L 2 means we are working with the L 2 norm, and the first index implies e.g. I\. 55 Chapter 3. Qualitative results =>: Because v\ G dli[u) and h[u] is Gateaux-differentiable, we have / 1 M - / 1 H > (tu-u,-ui), Vtu e F/,0, (3.21) i.e. because <9ii[u] is single valued. We know that v\ + v2 G d(h[u] + / 2 M ) , thus h[w] - h[u] + h[w) - h\u] >{w-u,v1+ v2) G VIfl. (3.22) Adding (3.21) and (3.22) we get: h[w] - I2[u] >{w- u,v2) Vwe V7,0. (3.23) Therefore, v2 G dl2[u). <=: By assumption dl\{u] is single valued and v\ G dl\[u}. Therefore, because v2 G dl2[u], (3.22) and (3.23) are satisfied. Hence vx + v2 G <9(/i[u] + I2[u\). • Thus, it follows that for fixed u G V>,0, Vw G VIfi: I[w] > I[u] - f{w- u)V • (S[u] + f) dO, (3.24) Jn i.e. v G dlhi[u] 4 « E —V • (S[u] + f). 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 p r o b l e m ' T h e o r e m 1 [Steady state problem] There exists a unique solution * s G Vj to (3.25), where ij) = \j>* _|_ U g ) a n d U s is foe minimiser of: inf I[v}. (3.26) vevt]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 results 3.1.4 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. • (3.27) Because pmin < p(c) < pmax, E is an strictly elliptic operator with bounded coefficients, thus for v G i7_1(ft), this problem has a unique solution z € H^(ft), (see [29], Th. 8.3 and Corollary 8.7 pg. 171), and we define E~l : if-1(ft) H-> HQ(£1) as this solution, i.e. £ _ 1[?J] = z, i.e. if p = constant, then i? - 1 = A - 1 . Note that because V-(S[u]+f) is the characterization of dlip.[u], this means that V-(S[u] + f) C L 2 ( f t ) , which is compactly embedded in iJ_ 1(ft). From (3.24) we can write: /(u;-u)V-'(S[u]+f)dn = [ (w-^EE^V-(S[u} + f) dft, JQ JQ = / (to - u)V • [pV(.E_1V • (S[u] + f))] dft, JQ integrating by parts we have: f(w- u)V • (S[u] + f) dft = - / V(io - u)pV(E-1V • (S[u] + f)) dft, JQ JQ because of the definition of E~l. 57 Chapter 3. Qualitative results We may note that since 0 < pmin < P < Pmax-, the associated seminorm defined by the inner product, < x,x > (v,w)vpQ= / pVv-Vwdfl, V u , w e V i o , Jn is equivalent to || • because of the zero boundary conditions. Note that this is just a weighted seminorm in HQ(Q) which we shall denote || • ||v"po-Therefore the subdifferential of I[u] with VPio seminorm and inner product is defined as: 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 ( £ - 1 V • (S[u] + f)) df2, Jn i.e. v G dlfy [u] => v £ E-lV • (S[u] + f). P ropos i t ion 5 If u e D ( d l i 2 ) then u G D(dlvpfi)-Proof Suppose it G D ( d l i 2 ) ; which implies that 3w G dlj?[u]. Thus I[ty] > I[u] - (w- u)v dtt \fw G V>)0 Jn =» /[to] > 7[u] - f ( w - ^EE^v dfl Vio G V/,0 Jo =*• J[io] > J[u] + / pV(w - u)V(E_1i;) dQ Vu; G VIFI Jn Because of the definition of the operator E[z] we have, I[w] > I[u] + / pV(w - u)Vz dtt Vto G Vifl. Jn 58 Chapter 3. Qualitative results Therefore z G dIvPi0[u] and because v G dl^[u] we can characterize z by, z G E^V • (S[u] + f). By definition, u G D(dl), provided that dl[u] ^ 0. Because 3z G dlvp0[u] this implies that u G D(dlvpfi)- • 3.1.5 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 function u G C([0, co); H) with u' € L°°(0,oo;H), such that, 1. u(Q) = u0, 2. u(t) e D(dl) for each t > 0, anal 3. -u'(t) 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 0 . We showed that dILi\u] ^ 0. Then, by Proposition 5 we have that if u0 G D(dIL2) then u0 G D(dIVf):0)-59 Chapter 3. Qualitative results Therefore, applying Theorem 2 to this initial condition, we have that there exists a unique function such that, 1. u(0) = u0, 2. u(t) e D(dIvPt0) for each t > 0, and 3. -u'(t) £ dIVp,0 for a.e. t > 0. This implies that -u'(t) £ E~lV • (S[u] + f), thus 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 rheo-logical properties. Behaviour of J ,^ Proposition 6 For a fixed gap-averaged speed, |V*|, and given £ {TY,I,TY,2\, then u £ C([0, oo); Vjfi) with u' £ L°°(0, oo; V7,0) -u'(t) £ E - 1 ? • (S[u] + f) -E[u'(t)] € V • (S[u] + f). (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 A1 fAH . / X J A'H2j(AH) 1 fAH ,i. . 0 = / ^(r)dr + f >- + r7'(r) dr, ' fAH A 1 M# 4' fAH J 1 /-A.H _ / r 2 _ , ( r ) d T + _ y T y ( T ) < J r , (3.30) where the ' denotes differentiation with respect to q. Thus, — fAHr7'(r) dr , A' = A[ ^ / V ' . (3.31) Jry Differentiating the constitutive law with respect to ry r(j) = Pool + Kjn + Ty, at fixed r, we find: <97 1 dry Poo + KUJ ,n—1 ' note that this is true for every r £ (ry,AH]. Differentiating the constitutive law with respect to r we have: fry _ 1 dr poo + Knin~l' Let / ( 7 ) = Moo + Knjn~1, then A' < — (3.32) ry x' + ^ < - + 4 (3-33) H Ty H (3.32) holds because r(7) > ry. Therefore X'<— '•• • (3-34) ry 61 Chapter 3. Qualitative results Proposition 7 For a fixed gap-averaged speed, | V * | , and given K* e [KI,K2]> then OK K tl K Proof Differentiating the constitutive law with respect to K, r(j) = Pool + K-yn+ Ty, at fixed r, we find:. ? j = 1! <n , 8K POO + Knjn 1 note that this is true for every r e (ry,AH]. Recall that, dj 1 Or poo + rcn-y™-1' As in Proposition 6, let f(j) = Moo + Kn'yn \ and note that Thus, T->r. K A' < — (3.36) K ^x' < - + ^f •• (3-37) K HK 62 Chapter 3. Qualitative results Proposition 8 For a fixed gap-averaged speed, | V * | , and given p^ € [p oo,iiMoo,2]j then d x < + (3.38) dpoo Moo Hp, Proof Differentiating the constitutive law with respect to n, T(i) = Pool + Kjn+ TY, at fixed T , we find: _ * L '= 2 < o, dpoo Moo + ^"7™ 1 note that this is true for every r 6 (ry, AH]. The proof follows in an analogous way as in Proposition (7), A' < — (3.39) Moo ^ x ' < J L + WL- •• (3-4°) Moo Moo Proposition 9 Let € H 2 , then (\A\ + \B\f-p(\A\p-2(A, A-B) + \B\P~2(B, B - A)) > C\A - B\2 (3.41) forC < 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\2 = (A, A - B) + (B, B - A), and because (A, A — B) > 0 and (B, B — A) > 0 we have: (A,A-B) = \A\P-2(A,A-B)\A\2-P <\A\P-2(A,A-B){\A\ + \B\)2-P (3.42) (B,B-A) = \B\P-2(B,B-A)\B\2-P<\B\P-2(B,B-A)(\A\ + \B\)2-P (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 b>a+d=a+ cosac > a + —, where a = rr — a. Thus, (|i4| + |B|)2-P(|i4|P-2(i4, A — B) + \B\P~2(B, B - A)) \A-B\ = (a + 6)2-p(ap-1 cos a + cos /?) > (a + 6)2-p(6p"1 cos /? - a p _ 1 cos /3) = (a + &)2"p(& - a)(p - l )^ - 2 cos /3 for b e [a, 6] = (a2-p + (2 - p)61-p&)(& - a)(p - l )^ - 2 cos/3 •> (2-p)62-p(6-a)(p-l)6p-2cos/3 (fr-a)(2-p)(p-l) 2 > JVc, where TY = 1/4(2 - p)(p - 1) < 1. Case c) /3 > TT/4, 7 < TT/2, Figure 3.2c. 65 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 — ^72 cos2 7 + sin 0:7 ~ cos a + sin 57 because 7 is small. Then, (|i4| + \B\)2'P(\A\P-2(A, A - B ) + \B\P~2{B, B - A)) \A-B\* (a + 6)2-p(&p-1 cos/3 - a?"11 cos a|) > c (a + b ) 2 - P / V ( b P - 1 7 + ( y 1 - a?-1)! cosa|; > -c sin 7 c Nb sin a b N > ~ 2 where ./V = 1/2. • Proposition 10 Let *2 <£ ^ (O), 6e too solutions for (2.63)-(2.65) for concentration fields c\ and c2 respectively, then j ( ^ j ^ f (V*i) - ^ | p ( V * 2 ) ) • V (*i " *2) d n * C||V(*i - * 2 ) | | £ 2 ( n ) , (3.44) /or a constant C which does not depend on $ior \&2. Proof 66 Chapter 3. Qualitative results Now, assume that exists f i + C fi such that V * i • V ( * i - * 2 ) > 0 and V * 2 • V ( * 2 - * i ) > 0 in f l + then, using Proposition 2 we have Therefore, x ( i y * i | ) v _ + x ( | v * 2 i ) v ^ _ ^ ^ c | v ( ^ _ ^ 2 ) | 2 _ ( 3 4 5 ) | V W i | | V W 2 | Without loss of generality assume that V * i • V ( * i —*2) < 0 in some region £l~ c fi, and just for notation let us write x i = x ( |V* i | ) and x 2 = x ( | V * 2 | ) thus, X l r V * i • V ( * i - * 2 ) + 7=^ 7-7 Vtf 2 • V ( * 2 - * i ) = . | v * i | 1 v , 1 | v * 2 | + ^ ^ . ( i v ^ r v f r i • v ( * i - * 2 ) + i v * 2 r v * 2 • v ( * 2 - *i)) + | V * 2 | - e V * 2 - V ( * 2 - * i ) ) f J X 2 1 X l 2 | V * 2 | 1 _ £ : 2 | V * i | ! - £ 1 X2 2 |V* 2 l - e | v * i | - £ v * i • v ( * i - * 2 ) + | v * 2 r e v * 2 • V ( * 2 - *i)) I V ^ I - V * ! • V(*! - *2)) (Ij^ pj -\^^) (3-46) for 0 < e < 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, , 6/Joo X " 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 * | ~ C 7 X m + 1 , thus, x ~ C | W | 1 / ( m + 1 ) c L—iv^i 1 /^" 1 " 1 ) - 1 , then X ~m + V X ' - (1 - e)-£- ~ CIV*! 1 /^ 1 ) - 1 (—!—• - (1 - e) which is greater than zero if (3.49) £ > 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 >Bmin-(l-e)Kmax, • (3.50) 68 Chapter 3. Qualitative results which is greater than zero if e>l-—. Ivmax Therefore, X |V#|i-e is an increasing monotone function for some range of e sufficiently close to 1. Using the assumption that V * R V ( * i - * 2 ) < 0 in 0~, we have that | V * 2 | > | V * i | in Or, then we have X L R V * i • V ( * i - * 2 ) + T ^ - T V * , • V ( * 2 - * i ) > |V*i | 1 v 1 " |V* 2 | X ( X l + (|V*i|- EV*i • V ( * ! - tt2) 2 V lV^ i l 1 " 6 |V*2|1 + i v * 2 r £ v * 2 - v ( * 2 - * i ) ) Using (3.5) in Proposition 2 and Proposition 9 we have, 1 / Xl , X2 (3.51) 2 V i v ^ i i 1 - ^ • IV*!!! 1- 6 v * i r £ v * i - v ( * i - * 2 ) + | V * 2 | - £ V * 2 • V ( * 2 - *!)) > c { | | ^ ± | | | | p | V ( * i - * 2)| 2 (3.52) Therefore, the result (3.44) holds. • 3.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 p2, while the other properties remain unchanged. We would like that our solution * i also varies continuously to *2, * i being the solution for p\ and * 2 the solution for p2. We present below a number of such continuity 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. Qualitative results 3.2.3 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(c2) u = ui, then: Jn ( * j w 2 f ( W 2 ) ) ' V ( U 1 " U 2 ) + V * 2 + V M l 1 _ 1 V * 2 , ) + f 2 ' V ( M 1 ~ U 2 ) d Q ~ °" ( 3 ' M ) Now, let v = U2 and u = u\, then: Jn ( X [ W ! | l ) ( W l ) ) ' V ( " 2 ~ " W l ) + + V w 2 l"IV*i | ) +fi • V ( « 2 - «i) dft > 0. (3.55) Note that, /^(|v*5 + v«i|)dn < /"^(|v*5 -v*n + |v*i|)dn (3.56), [ TF( |V*I + Vu2|)dft < / ^ ( | V ^ - V * t l + |V*2|)dft. (3.57) Jfi -n Jn n 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|) + (f2 - fi) • V(ui - ua) dfl > 0. (3.58) 70 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 /a(i^ <^>-w(^»)V(--»2)d^ /fl(W^ ">-2i^ fH-v(''-,")d0 + (X|v*2ll)(V^2) ~ X l v 7 1 l ' ) ( W l ) ) • V ( * 2 ~ ** } d Q - ( 3 ' 5 9 ) Similarly, / (f2 - fj) • V (m - ua) df i = / (f2 - fi) • V (*i - ¥ 2 ) dfl + [ (f2 - fO • - *I) dft Jn Jn Jn (3.60) Let us consider the following functional: 1 /.|V**+Vi;| x ( s i/2) ds 1/2 This functional is convex, lower semi-continuous, Gateaux differentiable and v ' | V * * + V u | v ' Using Proposition 5.4 in [22] we have: where Therefore, F ( « i ) - F ( u 2 ) >' (F'(u2), *i - *2) F(u2)-F(ui) > ( F ' ( « i ) , * 2 - * i ) (F'( U l), * 2 - * ! ) = jf ( ^ y f f v * ) • V ( * 2 - *i) df i . /n(wv*'-wS-v*-*!,d^0' ( 3 6 1 ) 71 Chapter 3. Qualitative results T h u s from (3.58) we have / (f 2 - f i ) • V ( * i - * 2 ) d f i + / (f 2 - f i ) • V ( t t $ - tfj) dfl + f 2 ^ | V ( * ^ - dfl Jn Jn Jn n + / X 1 ^ ( W 2 ) " l ^ ( W l ) ) - V ( ^ ~ * t ) d f t -Jn\ | V * i x ( | V * i | ) W i _ x ( | V ^ 2 l ) W 2 | . y ^ . ^ ) d o > Q . (3.62) U s i n g Cauchy-Schwarz inequality, A l & A 4 , the L i p c h i t z cont inui ty of * * w i t h respect to p, the facts that •vf lV7 \T/h and # i , # 2 € - f f 1 ^ ) . (3-62) becomes - P2\\LI > JQ v * i - * | v f f V * 2 ' ' V ( U l " U 2 ) d " - °- ( 3 ' 6 3 ) Therefore, using the result i n P ropos i t i on 10, we have | | P I - P 2 | | L * > C | | V ( * i - * 2 ) | l i a . (3.64) Yield stress continuity A s stated i n section §2.4.2 we know that x = x(l v"*|; H, Ty, K, n, Moo), thus we w i l l denote Xfc( |Vt f | ) = 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\ > Ty2, doing Taylor expansion for x i ( | V * | ) around 7y2 and neglecting 0((TYI — TYI)2), we have: x ( | V * i | , r y i ) = x ( | V * i | , r y 2 ) + ^ ; x ( | V * i | , r y ) ( 7 y i - r y 2 ) . (3.65) L e t t i n g v = u\ and u = u2 i n (3.53) we have: / (* 2 , ( jX*, 2 | ) V ^ a ) • V ( m - ua) + ^ ( | V M ^ + 7m | - | V * 2 | ) + f • V ( m - ua) dfl > 0. (3.66) Jn V |vw 2 | / w 72 Chapter 3. Qualitative results U s i n g the tr iangle inequality, the right hand side of (3.66) is less t han or equal to jf ( ^ ^ | ^ ^ p V ^ 2 ^ • V ( W l - u 2 ) - f - ^ ( | V * a - V * t H - | V * i | - | V * 2 | ) + f - V ( u i - u 2 ) d f l . (3.67) N o w let v = U2 and u = u\ i n (3.53)and use the tr iangle inequal i ty as i n (3.67), thus (^^p V*l) • V ( u 2 - U l ) + ^ ( | V * 5 - V * I | - r - | V * 2 | - | V * l | ) + f - V ( U 2 - « l ) dfl > 0. (3.68) A d d i n g (3.67) and (3.68) we get: + / ( ^ 2 - T V l ) ( | W i | _ | W 2 | ) d n Jn ti + [ ^ ^ ^ ( I V ^ i - V ^ I ) dfl>0. (3.69) Jn H U s i n g Cauchy-Schwarz inequal i ty and the fact that ry £ L°°(fl), the L i p c h i t z cont inui ty of * * w i t h respect to r y and the fact that \T/ € H1^), we can find upper bounds for the second and the t h i r d integrals on the left-hand side of (3.69). t (ZV2_ZiVi)(|V¥ 1 | _ | V g , 2 | ) dci < / " i T y a - T y i U V ^ i - V ^ l d n Jn H Jn < ||TYI - T y 2 | | L 2 | | * i - * 2 | | f f i < C||TVI - 7 v 2 | | i ~ ' (3-70) f ( T 2 2 ± Z n l | v * i - V * 5 | dn < C f | 7 y 2 - T y i | ( | T y 2 + T y 1 | ) dVl Jn H Jn < ||TVI - T y 2 | | i 2 | | T y i + r y 2 | | L 2 < C||7V1 - 7V2 | |£a < C||7Vl ~ 7 V 2 | | i ~ . (3-71) Subs t i tu t ing (3.65) into the first integral of the left-hand side of (3.69) we have: 73 Chapter 3. Qualitative results l i m - mlk2| |*i - *2|liji + C | | 7v i - m | | i 2 > [ ( M i n , - • - *,) df! > 0. (3.72) Jn V |vwi | | V W 2 | / 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: ^ ( i W 2 | ) ( v * 2 ) - 4 £ ^ ( v ^ | V * 9 | V * i | | V * I - * 5 | d f i + f I m - m l Jn dry |V( i t i -u 2 ) | dft. Using the result in Proposition 6, we get / Jn m , i - m l dry ;i - m l x( | v* i | ) + - |V(wi -u 2)| dft | V ( U 1 - U 2 ) | d f t<C / |7V,i Jn < C [ | 7 Y , i - 7 y , 2 | x ( | V * i | ) | V ( u i - U 2 ) | dft + C / I -TV.I -7y ) 2 | |V(u i -«2) | dft Jn Jn < C f | m - 7y , 2 |x( |V*i | ) |V(*i - * 2 ) | dft + C f |ry,i - TV, 2 |X( |V*I|) |V(*I - dft Jn Jn + C /" K , ! - r y , 2 | |V (* i -* 2 )| dft + C / | 7 y , i - 7 y , 2 | | V ( * I - ^ ) | d n (3.73) Jn Jn 74 Chapter 3. Qualitative results Using Holder's inequality and the fact that Ty^,min < 7Y,fc < 7Y,fc,max for k = 1 , 2 , recall that the rheological properties approach a constant value as £ —• 0 and £ —* Z, we have / i T V . l - T V . a l x d V ^ l D I V ^ - ^ l d f i < | | 7 V , l - 7 Y , 2 | | L = o | | X ( | V * l | ) | | L 2 | | V ( * l - * 2 ) | | L a 7 f i < C | | 7 y , i - 7 Y , 2 | | L » ( 3 . 7 4 ) because *fc e i 7 1 ( f t ) and x is a C°° function of |V*|. Similarly, / \TY,I - 7y , 2 | x ( |V*i | ) |V(* i - dft < ||TV,I - ry,2|||oo||x(|V*i|)||L2 < C\\TY+ - 7y , 2 | | £oo . J n ( 3 . 7 5 ) because the Lipchitz continuity of **.. Now, see that |V*J - *5| dft / Jn » ( l ' « . l ) ( m _ 2 a ^ l i ( v , 1 , |V* 2 | v ' |V*i | V* 2 | { 2 ) | V * i | < C | | 7 V i - T y 2 | | L i » < C | | 7 y i - 7 y 2 | | L ~ . ( 3 . 7 6 ) , ^ M X 2 ( | V * 2 | ) / V 7 , T , , X2(lV*i|) / V 7 , T , ,, ^ Assume that I ITVI — 7V2||L°° < i C l then |TVI - ry2||L°° > | | m - m i l l 0 Using this fact and combining ( 3 . 7 0 ) , ( 3 . 7 1 ) , ( 3 . 7 4 ) , ( 3 . 7 5 ) and ( 3 . 7 6 ) we have: l i m - T Y a l l i c c > jf ( X 2 ^ 2 l ) ( V * 2 ) - ^ ^ p ( V * ! ) ) • V(*i - *2) dft > 0 , using Proposition 1 0 we have the continuity result, U m - T r e l l i o o > C||V(*i - *2)||12. ( 3 . 7 7 ) 7 5 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, - « 2 | | i o o > C7||V(*i - * 2 ) | | i a , • (3-78) HAW - MOO,2||L~ > c | | v ( * i - * 2 ) | | £ a , (3.79) 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 = u2 in (3.80), then letting v = u2 and u = u\ and adding the resulting variational inequalities we have: J p 2V*2t • V(ui - u2) dfi - J piV*it • V(wi - u2) dfl + l(w< »^-i#(v*'»)V(«'--»df! + / 2^(|V*| - W?|) + (f2 - fi) • V(«i - ua) > 0. Jn H (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 dft [ p2V*2t • V(ui - u2) dfl- f pi W u • V(«i - u2) dfl Jn Jn = / P2(V*2t - V*u) • V(ui - u2) dft + / (p2 - pi)V*u • V(ui - u2) Jn Jn = I (P2 - Pi)V*u • V(*i - *2) dfl- I p 2(V* u - V*2t) • - tf2) dfl Jn Jn + f(P2- P i )W u • V(*^ - *J) dft - / p 2(V* l t - V* 2 t) • - *I) dft J Q Jn - / (P2 - Pi)V*« • V(* x - *2) dft - I f p 2^|V(*! - * 2)| 2 dft J O z Jn ot + f (P2 - Pi)Vtfit • V(tf$ - *I) dft - / p2(V*n - V* 2 t) • V(V2 - *I) dft (3.82) J S 2 Jn Combining (3.62) and (3.82) we have: ^(p2 - Pi)V* l t • V(*i - *2) dft - i j f p 2^|V(*i - * 2)| 2 dft + / (f2 - fi) • V(*i - tt2) dft + / 2^(.|V*5 - V*I|) dft Jn Jn w + y (p2 - p i ) w u • v(*^ - * D dft - y p2(v*it - w 2 i ) • v(*^ - * j ) dft > => / (P2 - Pi)V*« • V(*i - *2) dft + / (f2 - fi) • V(*i - *2) dft Jn J n + / 2^(|V*^ - WJ|) dft + / (ps - pi)V* l t • V(*5 - *I) dft Jn Jn - y p 2 ( V * i t - V * 2 0 - V ( ^ - * t ) d f t > X ^ | | V ( * i - * 2 ) | | | 2 (3.83) Because of condition, A.l p(c) e jL°°(ft). Having this in mind, the fact that <= ^(fl), ifr'k(t) e L°°(0, oo, if1^)) for fc = 1,2, the Lipschitz continuity of ** with respect to the 77 5 Chapter 3. Qualitative results physical parameters and using Holder's inequality, we have: ' / (p2 - Pl)V*U • V(*l - *2) dfl < | |p l -p2 | | i«» | |V(* l -*2)||£2 | |V* l t | |L2 Jn < C\\pi - P2\\L°° [(P2 - Pi)Wit • v(*5 - *i) dn < lbi-P2||i<»||v*iT||L2 < C\\pi -P2III0C • / (P2 - Pl)V(*l,t - *2,t) • V(*5 - *l) dfi < - P2|lioo||V(*l, t -Jn < C\\P! -P2WI00 [ ( f 2 - f i) • v (* i - * 2 ) dn < c u p ! - p 2 | | i - l i v ( * i - * 2 ) l l i 2 Jn < C | |p i — p 2 | | £ o o • /2^( |v*5-v*i | )dn < c||pi -P2 IU00 Jn -« Combining all this we have, 11 Pi " P 2 I U - > *^ | |V(* i - * 2 ) l ! | 2 (3.84) Using GronwalPs inequality in (3.84), we have the continuity result: | |p i-P2 | |L-T + | | V ( * 1 - * 2 ) | | i 2 ( 0 ) > ^ | | V ( * 1 - * 2 ) | l i 2 (3.85)-Yield stress continuity Using the variational inequality as above, for two different yield stress fields, ryt\, Ty,2, we have / p(W 2 t - Vtfit) • V(tf 1 - tf2) dil + [ p(V*2t - W u ) • V(*J - *3) dCl Jn Jn + .11*1 - ry2||L2 > Cjf (^^W) - ^ p ( V * i ) ) • V ( ? 1 - *2) d<7 > 0, (3.86) using the result in (3.77). Using the fact that ** is Lipschitz continuous with respect to TV and that *'(£) e H1, we can find an upper bound for the second integral in (3.86), thus 78 Chapter 3. Qualitative results ( p(V*2t - V*u) • V(*i - * 2 ) d f i + / p(V*2t - V*it) • V(*J - dO - " i l 1 1 ^ * 1 ~ * 2 ) i ' 2 l 2 + c l l r l y ~T2yll i2 (3-87) Combining the results in (3.77) and (3.87), we have the following differential inequality: l iny - T2YI\L2 > C^||V(tfi - * 2 ) | | | 2 (3.88) Using Gronwall's inequality in (3.88), we have the continuity result: l iny - r 2y|| L 2T + | |V(*i - *2)||x,»(0) > C||V(*i - *2)\\h (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, - K2\\L2T+ ||V(*i - #2)|M0) > C||V(*i - * 2 ) | | £ 2 (3.90) ||MOO,I-^OO,2||L2T + ||V(*I-*2)||L2(0)> C | | V ( * i - ^ H ^ . (3.91) 3.3 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. Qualitative results 3.3.1 End 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: Pin = TY,in = 7V(1), « i n = K ( 1 ) , " in = " ( 1 ) , Moo.tn = M o o ( l ) Pout=P(0), lY,out = T V ( 0 ) , K o u t = K ( 0 ) , n o „ t = n ( 0 ) , Moo.out = Moo(0) . The end condition ^in(<f),t) is the stream function that corresponds instantaneously to the steady parallel flow of a fluid with properties (pin,TY,in, Kim "in j Moo.in), 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: <9\I> 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, / ioo) , and therefore write: f1 <9\I>- f1 Q = %n(l)= _ ^ d 0 = / | V * | ( i 4 ; f f , 7 y i i n , K i f l ) n i n , / i o o , i n ) d .^ (3.92)' Jo Cl(p Jo 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 = Ain the solution of (3.92), and now define \I>jn as follows: r<t> Q\$ • r<t> *in(</>) = / -^rr(i(t>= |V*|(A i n ;H .Ty . in, Ki n ,n i n , /xoo.in) d0. (3-93) Jo ci<p Jo 8 0 Chapter 3. Qualitative results Similarly, *Out(0, t) is the stream function that corresponds instantaneously to the steady paral-lel flow through an infinitely long annulus of a fluid with properties (Pout,TY,out, Kout, nout, Moo,out), at flow rate Q ( t ) . We then define * O U i in exactly the same way as we have * j n , but using the rheological properties (ry)0Ut, K o u t , nout, Poo,out) • 3.3.2 Steady state problem In the classical setting, we characterise the steady problem as the solution * s G Vj of 0 = V-(S + f). (3.94). As discussed previously, this may be formulated more precisely as * s = ** + us, where us G VQ is the solution of the following minimisation: inf I[v], (3.95) v€V0 where I[u] is defined in (3.14). Note that if we choose ** = then evidently us = 0. Equally, from our consideration of the subgradient of I[u], we have that Vw G Vb: I[w] - I[ua] > - J(w - us)V • [S[u,] + f) dft. Multiplying (3.94) by \&s, integrating over ft, and using the divergence theorem: 0 = f * s S[uJ-«/ds+ / * S V - f - [ x ( | V * s | ) + 7 y / f r ] | V * s | dft, J a n J n (3.96) from which we see that I(uB) = - f L 4 ( | V * S | ) dft + / VsS{us} • v ds (3.97) J n J a n where i r\™s\* X ( s 1 / 2 ) L 4 ( | V * S | ) 1 r i v v s r v f s ^ l X(|V*S|) |V* S | - - JO ^ T T j ^ d s > 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 monotonically with |V*s| and is strictly positive for |V*g| > 0. 81 Chapter 3. Qualitative results Proof Taking the derivative of I/4(|V*S|) with respect to |V* S | , we have: L 4 = x'(|v*s|)|v*a| >0, with a minimum when |V* S | = 0. Therefore L4(|V*S|) is an increasing function and is strictly positive for |V* S | > 0. • 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 c2. These may for example correspond to different stages during a displacement. The two corresponding steady solutions are denoted ^s,i and *s,2-We 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, and tyout, are attained by and ^s,2 at each end of the annulus. Therefore, in homogenizing the boundary conditions, we may take the same ** for each concentration field, ci and c2. The consequence is that if Vs,i = * * + 115,1 and if ^s,2 = ** + uS,2, then ua>i is in the test space for «s i 2 , and vice versa. Lemma 2 Under assumptions AI-A4, I[u8,i,ci\ < I[uS,2,C2\ for any of the following situa-tions: 1. ry(ci) < ry(c2) a.e. (<£,£) G fl. 82 Chapter 3. Qualitative results 2. K(CI) < K(C2) a.e. (fag G ft. 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 with respect to any of the above. For example, let us assume that if ry (c\) < Ty(c2) then [X(|V*|) + TY/H](CI) < [x(|V*|) + TY/H](C2). Then we have: I[usA,ci] = inf I[v,5i] < inf I[v,c2] = I[usfi,c2\. veVo veVo This result then follows, provided that x(|V*|)+ry/i? is monotone in the directions indicated, for fixed IV*! Now A = x(|V*|) + ry/H is given implicitly by: 1 fAH • ' V * ' = ~A* J ^ ( R ) At fixed |V*| we differentiate with respect to any parameter q: -2A' fAH . . . . A'H2j(AH) 1 <~AH ° = " 7 5 - ^ 7 ( r ) d r + 1^-L + - J ^ r7'(r)dr, ^ / r 2 - 7 ( r ) d r + ^ / r 7 (r) dr, where the ' denotes differentiation with respect to q. The first integrand is positive, and therefore rAH sgn(A') = -sgn ^ r7'(r) dr^ . Differentiating the constitutive law, at fixed r, we find: dj 1 <o, dry Moo + « n 7 n _ 1 di 7 n <o, dn Moo + nnjn~l di 7 < 0. dpoo Moo + «;n 7 n _ 1 83 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 4(|V* 5 |)dO+ / *sS[u8] •»/ds Jn Jon 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 S ^ i d t = [L-^ + pcos[3dt = p(0)-p(L)+ fL pcosPdZ Jo Jo °^ Jo Therefore we have that: - /L 4(|V*s|,ci)dfi + gAp(c 1)<- / L 4 ( |W s | , c 2 ) dfi + QAp(c2), Jn Jn 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 — *s||x,2(n) decays as t —> oo. Afterwards, we compute numerically the solution for both systems and present the behaviour of ||* — \ & S | | L 2 ( O , ) - 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/S) and integrate by parts to get note that because = * s then (* - *s) e H^(fl). Now we multiply (2.63) by (\&a — \&) and integrate by parts, thus + Jn [ ^ ( | V * f l | - |V*|) + f • V ( * s - *) dft > 0. (3.100) Using the fact that x(|Vtt* + V"l) IV** + Vu| (|V** + V«|) is the Gateaux derivative of We can use Proposition 5.4 in [22] and have the following result: F ( * ) - F ( * s ) > <F '(* s ) ,*-tf s > where dfi. 85 Chapter 3. Qualitative results Thus, ( F , ( * B ) - f " ( * ) , * « - * > <0, because of Proposition 10. Therefore, adding (3.99) & (3.100) we have / pV^t • V(* s - tf) dfi - C||V(* - * s)| | | 2 > 0. Jn Because ^f- = 0, actually (3.101) states that -c||v(* - 9s)\\h > P~T|/n - * ) l 2 dn. Therefore, using GronwalPs inequality we have ||V(*-* s)|| £ a<||V(*-* s)|| i 2(0)e- K t ' (3.103) 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 if1(fi) for finite time. On the other hand, recall that numerically we work in finite dimensional subspaces of iJ1(Jl) -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 (3.101) (3.102) 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 1(f2), which is compactly embedded in W1,1+n(Q.). 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 ||* — * s | | L 2 ( Q ) 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 2 = 1, e = 0.4 and [3 = 0. Figure 3.4 b) shows the fixed concentration field, remember that we 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 Clearly this simulation agrees with the theoretical results stated in (3.103). 87 Chapter 3. Qualitative results Chapter 3. Qualitative results Chapter 4 S t a b i l i t y o f m u l t i l a y e r 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 multi-layer 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: — Both fluids are yielded at the interface. — Both fluids are 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 that fluid 2 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 single fluid along the annulus. We show that the flow is 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 parallel flow of 2 fluids. We consider the 4 types of base flow identified in §4.1. We show that linear instabilities can arise only if both fluids are 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, firstly for 2 Newtonian fluids and secondly for 2 power law fluids. 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 law fluids in 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 asymptot-ics of the stability problem. For long wavelengths it is possible to derive a semi-analytical expression for the eigenvalues. We also present an industrial application and make some comparisons with the results found in [59]. We find that interfacial instabilities are pre-dicted to occur only in the parameter space where the model in [59] predicts an unsteady displacement, i.e. a long finger on 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 perturba-tion 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 basic flows are parallel to the £-axis, hence * = *(</>). We assume that the annulus is occupied by two fluids: fluid 1 occupying <j> S [0,fa) and fluid 2 occupying <j> e {faA\- The interface is denoted <j> = fa. We denote the rheological parameters of fluid k as rky, rrik & Kk for k = 1,2 . Our solutions are steady parallel solutions of (2.83). Hence Sk = (Sk,4>,Sk^) = (5^,0), and (2.83) implies that: ^ 5 M = 0, (4.1) for the basic flow. By definition we have: _ dp pk cos (3 , . Sh*--Tz--sr-> (4-2) and with (2.89), this implies that the axial pressure gradient is continuous at the interface: |2 = 0. dp ~d~z 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 Sk^ is independent of <j> and £, (by definition, we look for solutions independent of £), following [59] we may write: SW = A, S2,(p = A - b , (4.3) 92 Chapter 4. Stability of multi layer Hows where b is the buoyancy parameter, given by: , P2~ Pi b = — cos 8. St* H (4.4) Note that typically, in a cementing scenario where fingering occurs, 6 is negative: the heavier 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, within fluid 1, and must be found as part of the base solution. In order to find A we use (2.85), with Q = 1 (for a fixed unit flow rate), d4> i = * r<t>i d$> r1 ( 1 ) = / d<p + Jo d<t> fc=i 4 dfa (4.5) fc=2. where, assuming the fluids are Herschel-Bulkley fluids, = sgn(A) { fc=i 0, H mi+2 /c111(mi+2) \A\ ~ TIX/H)MI+1 \A\ < T1IY/H, TI,Y/H \A\ + (4.6) 5* sgn(A - b) { fc=2 Hm2+2{\A-b\ - 7 i £ ) " i 2 + i K^{m2 + 2)\A-b\i mi + \A\ > ny/H, \A-b\< T2,Y/H, T2,Y/H-\A - b\ > T2,Y/H. (4.7) see [59]. Equation (4.5) is a nonlinear equation for A . It is straightforward to show that-the 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 two fluids. 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 the first derivative (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 basic flow depends on 9 parameters, rytk, «fc, mis, f° r k = 1, 2, b, e and </>j. The buoyancy itself depends on 4 parameters, pk, St* and /?. Although these do not appear in the base flow description, 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 parallel flow of two generalised Newtonian fluids with rheological parameters, Ty,\ = 2, «i = 2, m\ = 2, ry,2 = 1, « 2 = 1, wi2 = 1, 6 = —9.5, the eccentricity is 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, both fluids are fully yielded, the pressure gradient and buoyancy forces are large enough to over come both Ty\ and Ty2- The stream function is continuous and a discontinuity on the axial velocity is observed at the interface. Figure 4.2 shows the basic parallel flow of a couple of fluids with rheological parameters, ry,i = 0.5, Ki = 1, mi = 2, Ty,2 = 0.6, K2 = 1, m,2 = 1, b = —9.5, the eccentricity is 0.5 and fa = 0.5. The inclination of the well is (3 = TT/10. The pressure gradient and buoyancy forces are not sufficiently large to fully over come Ty2 along the narrow part of the annulus and a small mud 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, m2 = 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, m2 = 1, 6 = 0,the eccentricity is 0.2 and fa = 0.5. For this case we have a fully vertical 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 i 2 = 1, «2 =T, pi = 1, m 2 = 1, e = 0.4, 4>i = 0.5 and /3 = TT'/IO. I1 Figure 4.2: Multilayer base flow, both fluids yielded at the interface, mud channel in the narrow-side of the annulus. Rhelogical and physical parameters are: ryti = 0.5, K\ = 1, p\ = 2, m i = 2, 7Y I 2 = 0.6, «2 = 1, P2 = 1, m 2 = 1, e = 0.5, ^ = 0.5 and [3 = TT/10. 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 2 = 1, P2 — 1, m2 = 1, e = 0.8, & = 0.5 and B = TT/10. Figure 4.4: Multilayer base flow, fluid 1 unyielded at the interface. Rhelogical and physical parameters are: Tyti = 7, K\ = 7, p\ = 1, m\ = 2, T y ) 2 = 0.2, K2 = 1, p2 = 1, m2 = 1, e = 0.2, <?!>i = 0.5 and /3 = o'. 96 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 which fluid 2 is static, either partly or fully, i.e. there is a static mud channel. Let us consider a basic multi-layer flow of two fluids with 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 in fluid 2 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 dfa (4.9) fc=i io d<t>\ i.e. all the unit flow rate must flow through fluid 1, <f> G (0,fa\. Evidently as fa —> 0, if the mud remains static, then A(fa) —> oo since the unit flow rate is forced through an increasingly narrow azimuthal gap. It follows that for some fa, 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,min, and the maximal azimuthal width of mud channel is therefore 1 - fa min- Let us now explore the parametric dependency of (4.8) and (4.10). 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 d* k=l n,Y\MI \n,Y ( A H _ \ m i + 1 (AH + 1 \Ti,Y mi + 1 „ , ,'AH\ mi + 2 Jl,Y J (4.11) which at each <f> depends only on e, mi, Ajr\y and B\ = T\y/«i- Therefore, we see that (4.9) implies A(fa) T\,Y /(e,Bi,mi). Consequently, dividing (4.10) by r^y, we have H(fa^min) f(e,Bi,mi) n,Y T2,Y i.e. depends on (e,B\, mi, At, Ay), where Ah = ^ Av T"2,y (4.12) (4.13) T\,Y n,Y are the buoyancy ratio and yield stress ratio, respectively. The condition that (4.8) holds is now: Ay > (1 - e) /(e,Bi,mi) n,y (4.14) 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 <&,min- As expected, if we 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 ry;2 > ry,i, i.e. Ay > 1, the pressure gradient necessary to overcome the yield stress of both fluids at the interface will not suffice to overcome the yield stress of fluid 2 in the narrow side of the annulus. Having as a result a decrease of the domain where fa^m = 1- The same phenomena will occur if either n2 > «i and p2 > Pi-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 2 = 0.4, p2 = 1, rn2 = 1 and 3 = 0. The zero level curve corresponds,to the equality of, condition (4.14), i.e. < )^mjn = 1, the rest of the level curves correspond to different values of &,min when condition (4.14) is fulfilled. 6 <1 • ' Av-(1-e)|l(e,B],m])-ABl=0 Av-(l-e)|l{e,BI1mI)-Aj>0 No static tayer Static layer 0.2 0.4 0.6 Figure 4.6: Contour plot of 4>i,min, 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^min = 1, the rest of the level curves correspond to different values of r^.min when condition (4.14) is fulfilled. 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 <foim , 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, K2 = 0.4, P2 = 1.1, m2 = 1 and 8 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. <?!>iimin = 1, the rest of the level curves correspond to different values of <pi}min when condition (4.14) is fulfilled. Figure 4.8: Contour plot of 4>i<min, 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, p2 = 1, m 2 — 1 and 8 = 0. The zero level curve corresponds to the equality of condition (4.14), i.e. Si m i n = 1 the rest of the level curves correspond to different values of ^ ^ m i n when condition (4'.14) is fulfilled. 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 0.8 1 Ay Figure 4.9: Contour plot of falTnin, 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, m2 = 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 of fa}min when condition (4.14) is fulfilled. 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 to find \&. 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,mm, and for certain combinations of parameters. We have shown that fatmm depends only on 6 parameters, ry ,^ for k = 1, 2, «i, mi, b and e. Therefore, provided that fa > fa,nun 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-layer flows in §4.1. We denote the transient 101 Chapter 4. Stability of multi layer flows interface by, F = <p — fa(£,t) = 0. Our system of equations can be written as: V • [Sfc + pfcV*fc,t] = 0, mflk, (4.15) ~d!T ~dqV~dJ = ' ^4> = fa(gt) (4.16) *i(0,^,t) = 0, • (4.17) *2(U,t) = 1, (4-18), %lfc= 2 = 0, at <£ = <&(£,£), (4.19) (Sfc + PfeV*fc,t) • n|JL2 = -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 = (Sfci(^ o,0), with interface position fa = fap, and recall that Sk,<j>,o is constant in each domain. For <5 < 1 we assume the expansions: *fc = ^k,o + 5^k,i+S2^>k,2 + - , Sfc = (£fc,,/,,o,0) + JSfcj + <52Sfc)2 + ... fa = fa,o + Sh + ..., substitute into (4.15-4.20) and expand in powers of 5. V* f c = (*fc,o,* + ^ M , * + -.^fc,U + -) (4-21) |V*| = [(*fc,o,* + <J*M,* + -) 2 + (***,U + -) 2]' ~ I*fc,o,*l+^gn(*fc,o,^ *fc,i,^ +0(^ ). (4.22) When the fluid is yielded we find: Sfc,i = [xfc(|*fc,o,*l)^,i,*. (Xfc(l**,o,*|) + ! W 1 ) T | ^ M T ] ) (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 ^ , Kk, mjt, for fc = 1, 2, b, e and 102 Chapter 4. Stability of multi layer Hows <f>i. Now the buoyancy parameter b is defined in terms of the two densities, pk, the inclination [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 the flow of a single fluid along the annulus, and assume that the base flow of the single fluid is everywhere yielded, although the arguments following will also apply if the base flow is 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*!)t = V • x'(l*o,*|)*i,*, X(l*0,*|) + T7 (4.24) with boundary conditions: *x(0) = 0 = *i(l). We suppose a normal mode solution for * i : tf i ~ /(</>)ei(c*-st). Substituting \I>i into (4.24) we have: isp(D2 - a2)} = D[DfX'(\*0M ~ a-2X(f*0,*|)+ jy (4.25) / 103 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[Jx+a2h] = J2 + a2I2, (4.26) where -fi = /% l / l 2 d* (4-27) k _ / V * ( l * ¥ V * d * (4.28, JO l*O,0| Ji = / p|£/|2d</, (4.29) Jo h = I \ D f \ 2 X ' ^ (4-30) Jo 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 S j , the upper and lower analytical bounds are defined by: — maxx' < Si < —minx', (4-31) The argument of x' is |*o,<^ |(^ ) a n d the max and niin are taken over <f> £ [0,1]. It has been 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 = 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. 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(X + Ty/H) min|^| - 1 - max|*'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 is V • [Sfc,i + AbVtyfc,i„i] = 0, 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) * 2 , i ( U , i ) = 0, (4.36) 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,*)|fcI2 = 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) n, about (j) = 4>ifi w e have: [Sfc + PfcV*fc,t] • n ~ Sk,o,(j> + Sh^Sk,o,4> + 5Sk,i,(i, +Pfc[*fc,o,tft + Sh^kfiMt + s^k,i,4>t\ at <j> = fop (I,tan/3sin7r0j) • n ~ 1 — Sh^ tan 8sinnfafi 106 / 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)lfc=i = 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 £ [ X 2 ( l * w l ) I ? / 2 ] - a 2 isPl(D2 -c?)h isp2{D2 - o?)f2 0, 0, and at H sh0 + a(fk + */c,o,0^ o) = 0, (fto*o,M + /fc)li = °> [(Xk(\Vo,k,<j>\)-ispk)Dfk]\lJl = iabho tan (3 sin ir<f>ifi (4.39) (4.40) (4.41) (4.42) (4.43) (4.44) (4.45) 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> > fa we have 32,0 < H(4>iY By choosing 5 sufficiently small we may ensure that T2.Y |S2,o + 5S2,i| < 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: D[xi(|*i,o,*IWi]-a2 xi ( l*wl) + ! ^ r i fi, = isPl{Dz-al)h l*i,o,*l > 6 ( 0 , ^ , 0 ) (4.46) /i(0) = 0, ' (4.47) /i(&,o) ' = 0, (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 5fayt\(£>,t) is the linear perturbation of the yield surface. Note that fa < fay < 1 and fa < fa,min-As 4> —> fatY fluid 2 becomes unyielded across the annular gap and static, i.e. 1X7*21-* 0 as <j>-^fatY.. Expanding each component to 0(5) and linearising onto fayfl, we have that: * 2 , u + <A2,y,i(^i)*2,o,^ = 0, at <f> = fayfi, (4.49) * 2 , u = °' at (A = </>2,y,o. (4.50) Equation (4.50) leads to the boundary condition, / 2 = 0 at <j> = 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 stress fluids, [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? ^ r ^ - l / 2 , = isp2(D'-a*)h l*2,0,*l 4> G (fafi, fayo) (4.51) Hfayfi) = 0, (4.52) i.e. only the domain of the fluid 2 problem is changed. The other situation to consider is that in which fluid 1 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 basic flow. It follows that |Si | < Tiy/H at the interface, for the base flow, 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: i V111M/1 n^,i) -t-U[xi(|*W|)I>/i]-a2 Xi(l*i,o,*l + • R h , = ispi{D2-o?)h l*i,o,*l 0€(O,0i,y,O) (4.53) with boundary conditions (4.41) and /I(0I,Y,O) = O. (4-54) 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 base flows essentially reduces to the stability of one or more independent single fluid layers. 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 Newtonian fluids in a concentric annulus (H(<p) = 1), we find: . 6 «i b 1 + (1-</>i,o)~— 0»,OQ— *o^ = • % * W = * t ' V ' (4"55) 4>i,0 + (1 - <Pi,0)— <Pi,0 + (1 - <Pi,0) — K2 K-2 110 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 M - V i s ) C D 2 - a 2 ) / i (3«2 - ip2s)(D2 - o?)h 0, 0, <f> e (0,fa,o) <i> e (fa,o, 1) We note that if (ink — ipks) = 0, then sj < 0 and these modes are stable. Therefore, for instability we may assume that (3«fc — ipks) ^ 0, and interestingly observe that the eigenvalue 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 (4.41) & (4.42) are: fi{fa) =Ai sinha0, f2(fa) = A 2sinha(0-1), (4.56) for constants A\ & A2. We may assume that ho ^ 0, write A\ = hrjCit ^42 = hoC2 for constants Ci & C2, and divide (4.43)-(4.45) by hQ. From (4.43) & (4.44) we find 1 Ci C2 = 1 sinh a fa o 1 + a e fafi + (l-fafi) sinha(^ i,o - 1) « l «2 6 «2 3/t2 0i,O a e fafi + (1 - & , o ) «2 Finally, substituting into (4.45), we find the following quadratic relation for s: i tan B sin nfa o 3/t2 •cotha^n ac 1 + 0i,O + (1 - fa,o) • cotha(fafi — 1) K2 b K l • Pi Z S — K 2 3 K 2 K l « 2 3K 2 02,0 a £ 0i,O + (1 - 0i,O) K2 1 — is P2_ 3 « 2 (4.57) 111 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: K l = l,e = 0.0, fa = 0.55 and f3 = 0. a) px = 1, P2 = 0.75. b) pi = 1, p2 = 1.25. 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, the flow may be unstable if the heavier fluid is on top of the lighter one, i.e. p\ > p2. If we have positive buoyancy, i.e. p\ < p2, then the inclination has stabilizing effects. 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 each fluid layer. 4.5.2 Power-Law fluids in a concentric annulus Because the base velocity is constant in each fluid layer, we find: \A\mi \A — b\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,*l2|) = \A - b\. (4.59) Using (4.59) we can find an explicit expression for x'k m terms of A , , , . T ., <1(mi + 2) .. K2™2(m2 + 2) and Xi(l*o,*,il) = «r(mi + 2) X2(|*0A2l) = < 2(^ 2 + 2) |*o,*,i| l^l"11-1 l*0A2l \A-b\^-i • y • > 115 Chapter 4. Stability of multi layer flows Substituting (4.60) & (4.61), into (4.39) & (4.40), we get DVi_q2 /^*>-wry\ = 0, ,e(0,,(0) (4,2) V 1 21+)-Mmi-w D h-a — m = 0, 0G(<Pi,o,l) (4-63) with boundary and jump conditions: /i(0) = 0 (4.64) /2(1) = 0 (4.65) (4.66) • L * • , / K m i (m 1 + 2) \ /K 2 n 2 (m 2 + 2) \ (4.67) and the kinematic condition — + H / l + o s g n ( a ( / i /iosg (A) j T O l' ' , 0 , ) = 0 at 0 = &,0 (4.68) Solving (4.62) & (4.63) with (4.64) & (4.65) respectively we have: /1(0) = Ci/iosinhAir\ /2(0) = C2/i0sinhA2(</>-l), (4.69) for constants C\ & C 2 and A2 = ^ f ^ ( m i + ^ - t l ^ - V t S A 2 = a 2 / / t 2 " 2 ( ^ 2 + 2 ) - i ] A - 6 r - y 2 S 116 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 (s sg*(A)\A\mi' s i n h A i ^ j \ae KMI ( m i +2))' In a similar way, using the first jump condition, (4.66), and the definition of C\, we can find an expression for C2, C2 = -r s sgn(A-b)\A-b\m* sinhA20j \ae K2n2{m2 + 2) Using the definitions of C\ & C2 and (4.67) we find a non-linear equation for s, iabho tan 8 sin IT fa o «7 l l (mi+2) m i l ^ l ™ ! - 1 K 2 n a ( m 2 + 2) m2 A|m2-1 • pis P2S Ai tanh Ai<; A2 tanh A2<j s sgn(A)\A\* ae KMI(MI + 2) s_ _ sgn(A- b)\A -b\m* c7e~ K%2{m2 + 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, K2 = 1, e = 0.0 and 8 = 0, mi = 1.5, m2 = 2. a) fa = 0.5. b) fa = 0.6. For b < 0, p\ = 1 and p2 = 0.75, for b = 0, p\ = 1 and p2 = 1 and for b > 0, p\ = 1 and p2 = 1-25. 117 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, m2 = 2. For b < 0, p\ = 1 and P2 = 0.75, for b = 0, p\ = 1 and p2 = 1 and for 6 > 0, pi = 1 and p2 = 1.25. 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 as for the 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 results for the eigenvalue problem, D [ x i ( l * w I W i ] - « 2 £[x2(l*2,o>l)Af2]-a2 isPl(D2 - a 2 ) h isp2(D' - c O / 2 0, 0, and at Pi,Q-H sh0 + a(fk + ^k,o,4,ho) = 0, (Ao*o,M + Jfc)li = 0, k=2 [(Xk(\^o,k,<t>\) ~ ispk)Dfk] k = l = iabh0 tan j3sinirfa,0 (4.71) (4.72) (4.73) (4.74) (4.75) (4.76) (4.77) 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: Tn{y) = cos(n arccos(j/)). • As solutions of the following Sturm-Liouville problem, ( v^Ar . ( y ) ) + - ^=r B ( y ) -o . dy Vv " dy-"""J • As a recurrence relation, T0(y) = l, Ti(y) = y, ... ,Tn+1(y) = 2yTn(y) - Tn_i(«). 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) = J2anTn{y), and evaluate the Chebyshev polynomials at the extrema of the iV-th Chebyshev polyno-mial, 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) 2 -(-fa + <Pv)y + (4>i + 4>Y) for 0 6 [0,(pi], 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 = cBa. 120 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 f2{4>) matrices. In the case of boundary conditions, we implement them only on the first and 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 1 1 1 1 1 -0.21 1 ' ' 1 1 "0 10 20 30 40 50 0 10 20 30 40 50 a a (a) (b) Figure 4.18: Maximum imaginary part of s vs. a. Inclined annulus. Rheological and physical parameters are: «i = 1.5, K2 = 1, m\ = 1, m2 — 1, b = —1 and 3 = 0. a) Exact solution, b) Numerical solution. 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 1.635732937987910E-10 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.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 the flow in a concentric annulus. As expected, because both fluids are yielded at the interface, yield stress does not appear to have a large influence on the stability of the flow, and nor does interface position. Only viscous and buoyant forces govern the unstable behaviour. The results are very similar to those of the Newtonian fluids in Figures (4.11)-(4.14). Figure 4.20a) shows similar behaviour as for Newtonian fluids, compare with Figure (4.15). As above, both fluids are 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 lighter fluid 2 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 Tyti = 0.7. a) fa = 0.5. b) fa = 0.6. For b < 0, p\ = 1 and p2 = 0.75, for b = 0, p\ = 1 and pi = 1 and for b > 0, p\ = 1 and p2 = 1.25. 123 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, re2 = 1, e = 0.0, fa = 0.5 and 8 = 0', mi = 1, m2 = 1, ryj = 0.7 and ry]2 = 0.9. b) «i = 1.5, K2 = 1, e = 0.0, & = 0.5, /? = ir/8, mi = 1, m2 = 1, ry,i = 0.7 and ryi2'= 0.9., mi = 2, m2 = 1.5. For b < 0, p\ = 1 and p2 = 0.75, for 6 = 0, p\ = 1 and p2 = 1 and for b > 0, pi = 1 and p2 = 1.25. 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 desta-bilizing 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 2 = 1, e = 0.0 and 3 — 0, mi = 1, m2 = 1, 7y , i = 0.9 and 7 y 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 p2 = 0.75, for b = 0, pi = 1 and p2 = 1 and for b > 0, pi = 1 and p2 = 1.25. 125 Chapter 4. Stability of multi layer flows -6 .§i 1 1 1 -0.5 0 0.5 1 I i (a) 1| r- , , , 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, K2 = 1, mi = 1, m 2 = 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 f2, dashed line shows the imaginary part of fi and f2 and dashed thick line shows the magnitude of fi and f2. c) Eigenfunctions of an interior mode, solid line shows the real part of fi and f2, dashed line shows the imaginary part of fi and f2 and dashed thick line shows the magnitude of fi and f2. 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, K2 = 1, mi = 1, m2 = 1, Tyti = 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\ + a2S2 + ••• 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) ) Sl + e ( *0,1,* — (<t>ifi - 1) - *0,2,tf<fo,0 ) = 0 K2 I \ K2 Thus, e ^*o,i,* (^l " &,o) + *o,2,*&,o) S l = v r. Kl \ '- (4-78) 4>i,0 H (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, n2, pi, P2, St* and 4>i$- Now, we continue to the 0(a) problem. O(a) problem Collecting 0(a) terms, we have the following equation for s2, % \ • 2 ( P2 1 Pi 1 e\<pifl-l 4>iflJ \3K2 e(cj>ifi -1) 3/t2e<fo,0 \3K2<Pi,0 . 3«2 (<pi,0 - 1) 4- «~—tanpsin7r</>jo, 6K2 Therefore, P2 , S2 = ISi I Si I 0iO 1 \3K 2 ~ ^,0)^,0 tan/?sin7r ,^o *i,o + £(l-&,o) ' ( 7 9 ) Clearly, s2 is pure imaginary and in the presence of negative buoyancy, i.e. fluid 1 heavier and 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* - *o,2,*) + (*o,2,^ - ^.l,*)) > 0, (4.80) \Pl «2 / the inequality is reversed if s\ < 0. This will lead us to 2 cases, 1.) Case 1, thus (*o,2,* " *o,i,*) > 0, b < (4.81) «2 P2 2.) Case 2, thus (*0,2,* - *0,1,*) < 0, K l - «2 > jj ^ > (4.82) «2 P2 If S i > 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 the flow. It is worth to mention that this phenomena is not comparable to Kelvin-Helmholtz instabilities, see e.g. [18]. This classical instability involves two inviscid fluids one 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 the fluids is 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 parallel flow of two viscous fluids in 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 the flow is 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 the flow if the heavier fluid is 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 s2 with respect to p2 and /t2. We have 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, (*o,2,</> — *o,i,4>) > 0, from (4.55) we have, « 2 130 Chapter 4. Stability of multi layer flows and for instability to arise, from Case 1 we have, 2 ! < £ = 1. K2 Pi Therefore the flow will always be stable. \ \\ • \ c c \ o \ \ V ^ " / V? • / / / / / \\'8 / \ / Figure 4.23: Contour plot of s2. Rheological and physical parameters are: K\ = 1, p\ = 1, e = 0.0, fa = 0.5 and 8 = 0. 4.7.2 Yie ld 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 ~ as i + a s2 + ... h ~ + och\ + a2h2 + 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 and at 6 = D(x'kDfkfl) = 0 (4.83) /fc,o = 0 0 = 0,1 (4.84) H —jsiho + (fkfl + ^kfirfho) = 0 (/fc,o + *fc,o>o)li = 0 x'kDfkfilt = 0. .(4.85) (4.86) (4.87) Solving (4.83) we have: where /!,O(0) = Jo Xi /2,O(0) = /2,O(0i)+ / ^d4> J<t>t Xi C\ /2,o(0i) = / —r d0 + ^o*fe,o>l2-Applying the boundary condition at 0 = 1, we can find an expression for Ci, (4.88) (4.89) Ci /i(0O + /2(l) where 1 [<>> 1 L(0) = / - d 0 J2(0) = / - ^ 0 . JO Xl J(pi A 2 (4.90) (4.91) 132 Chapter 4. Stability of multi layer Hows From (4.85) we get clearly, si G R. 0(a) problem D(X'kDfk,i) = PkisiD2fkfi fk,i = 0 0,1 (4.93) (4.94) and at 4> = fa U — (s2ho + sihi) + (fk,i+Vk,o,<l>hi) = 0 (4-95) (/jt,i + **,o,^i)li = 0 (4-96) (x'kDfk,i-isiPkDfk,o)\i = ibfro tan/?sin71-^0• (4.97) Integrating (4.93) once, we have: x'kDfk,i-PkisiDfkfi = C2,k. .(4.98) using (4.97) we get, C2,2 = C2,i + ibho tan dsinn fao Integrating (4.98) and using (4.96) we can find the solutions to (4.93), 133 Chapter 4. Stability of multi layer flows /2,i(0) = h M - h i M l k M + j + i s i K ^ ) d ^ (4'10°) where C\ is defined in (4.90). Using the boundary condition at 4> = 1 we can find an expression for C 2 , i C 2 , i = i j / J . \ i T t -\ \ -C\isi T ) , \ , T m -t6feotan^sxn7r0i,o r , v (4.101) h(<t>i) + h(l) hWi) + h(l) Ii{4>i) +12{±) where I\ & I2 are defined in (4.91) and Using the kinematic condition (4.95), (4.99), (4.90) and (4.101) we have an expression for s 2, g(&,0) . , + . . . , Ji(^)/2(1) , ,„ A , T , Jx{fa)I2{l) - J2{l)h{4>i) , . --—s2 = , 6 t a n / W f r o f i W + f a ( 1 ) + M l A * w ( / l ( ^ ) + I2(l))2 ' ( 4 ' 1 0 3 ) clearly, s 2 = s2,/. Note the first term of the right hand side of (4.103), as in previous cases, if b < 0, i.e. heavier fluid 1 is on top of lighter fluid 2, the buoyancy has destabilizing effects. 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 11 parameters on (4.103). Fortunately, from §4.6.2 we know 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 the flow when we are in the case of a mildly eccentric vertical annulus. Figure 4.24 shows the level curves of the maximum imaginary part of s2 for rheological para-meters, pi = 1, pi = 1, m\ = 1, m2 = 1, Ty}\ = 1, Tyt2 = 1. The neutral stability level curve, Ki = K2 implies that if fci > k2, i.e. the displacing fluid is less viscous than the displaced one, the flow will be stable. In Figure 4.25 a) and Figure 4.25 b) we have decreased and increased, respectively, the buoyancy parameter. When the displacing fluid, fluid 1, is denser than the displaced fluid, fluid 2, 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 displaced fluid is 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 than fluid 2, resulting in a shift of the level curve for n\ = K 2 to the stable region. The opposite effect happens when mi = 1 and m2 = 1.2, Figure 4.26 b). Now the level curve for K\ = K 2 has been shifted to the unstable region, clear effect of fluid 2 being more shear-thinning than fluid 1. Figure 4.27 shows the effects of yield stress on the stability of the flow, for Figure 4.27 a) fluid 1 has a higher yield stress than fluid 2, TY,I = 1.5 and 7 y i 2 = 1. Because of this, a higher pressure gradient is required to overcome 7 ^ 1 for both fluids to be yielded at the interface. This could be seen as fluid 1 being "more viscous", thus the level curve K\ = K 2 is shifted to the unstable region. If now 7 y ) 2 = 1-5 and Tyt\ = 1, Figure 4.27 b), we see that fluid 2 is "more viscous" and the level curve «i-= K 2 is shifted to the stable region. To the knowledge of the author, this is a phenomena never reported in the literature. For this kind of flow, even in the absence of inertia, the flow can be unstable. All rheological parameters 135 Chapter 4. Stability of multi layer Hows play a role on the s tabi l i ty of the flow, something not seen i n Newton ian and power law fluids. For y ie ld stress fluids, even wi thout the presence of buoyant forces the flow can be unstable. Fo l lowing the results found i n [59] we can make some comparison of our contour plots and unstable regions for the cases when the displacement is unstable. F igu re 4.24: Contour plot of s2, « i vs. K2. Rheological and physical parameters are: p\ = 1, m i = 1, m 2 = 1, 7 V , i = 1, 7v ,2 = 1, P2 = 1, e = 0.2, fc = 0.5 and (3 = 0. 3| 136 Chapter 4. Stability of multi layer Sows Figure 4.25: Contour plot of s2, «i vs. K2. Rheological and physical parameters are: mi = 1, 77i2 = 1, TY,I = 1, Ty,2 = 1, e = 0.2, fa = 0.5 and 8 = 0. a)pi = 1.1 and p2 = 1. b)pi = 1 and P2 = l . l K1 (a) (b) Figure 4.26: Contour plot of s2, KI VS. K2. Rheological and physical parameters are: p\ = 1, TYi = 1, ryi2 = 1, p2 = 1.1, e = 0.2, fa = 0.5 and /3 = 0. a) mi = 1.2 and m 2 = 1. b) mi = 1 and m2 = 1.2, 137 Chapter 4. Stability of multi layer Hows 2.5 1.5 1 r 1 ? / u 0.01' J / S 1.5 2 2.5 (a) (b) Figure 4.27: Contour plot of s2, KI VS. n2, buoyancy parameter equal zero. Rheological and physical parameters are: K\ = 1, p\ = 1, K2 = 1, p2 = 1, e = 0.2, fa = 0.5 and 8 = 0. a) 7y ;i = 1.5 and Ty,2 = 1. b) Ty,i = 1 and 7y}2 = 1.5. 4.7.3 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 n2, i.e. viscosity of fluid 2 increasing, an increase of the domain of 138 Chapter 4. Stability of multi layer Hows unsteadiness and ins tabi l i ty is observed. In F igu re 4.29 we show the level curves of the m a x i m u m imaginary part of s2 as a funct ion fo TY,2 and e. W e have the same fixed rheological parameters as i n F igure 4.28 w i t h fa = 0.5 and 8 = 0. W e capture the same behavior as i n F igure 4.28, as fluid 2 becomes more viscous, the unstable doma in increases. Surprisingly, our unstable domain occurs on ly as a subset of the unsteady and fingering domain i n F igure 4.28. T h i s means that i f a long finger develops du r ing the displacement, for certain parameter ranges the paral lel flow may be unstable. Therefore, apart from having an unsteady displacement, m i x i n g may occur at the interface of the finger due to the ins tabi l i ty that we s tudy here. It is not clear why there is apparent ly no interfacial ins tab i l i ty for parameter regimes for wh ich the annular displacement would be steady. No te that for such parameters, a base paral lel flow is a perfectly well-defined solut ion of the two fluid problem. 139 Chapter 4. Stability of multi layer flows Figure 4.28: Example unsteady displacement and local instability/viscous fingering regimes: effects of changing ry,2 and e for different K2- (a) K,2 = 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). 140 Chapter 4. Stability of multi layer flows 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 (a) • \ • \ , \ \ \ J.0I s s •\ r ^ U i I 1 l_ 5.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 e (b) (c) Figure 4.29: Maximum imaginary part of s vs. a. Rheological and physical parameters are: KI = 1, 7V,i = 1, mi = 1, 7 7 i 2 = 2, 6 = —0.5, '= 0.5 and (3 = 0. a) K2 = 0.5. b) n2 = 1. c) K2 = 1.5. 4.8 Spatial stability problem Finally, we present a collection of results concerning linear spatial instabilities. Although in-complete, 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, for the single fluid problem, a spatially developing instability: s > 0 and a = aij + ia/. In this case, (4.25) constitutes an eigenvalue problem for a2 & C. We find: -(a2)R = a2j-a2R (a 2 ) / = lajaji Solving for aR h a/, we have: ai = ± aR = ± It appears that any non-trivial complex eigenvalue a2 and eigenfunction / must admit spatially growing instabilities. 4.8.2 Two fluids spatial instability In an analogous way as in the previous section, let us assume that: s > 0 and a = aR + ia/. We assume small s and study the asymptotic behavior of the problem. Thus, let /fc a h substituting this expansion into (4.39)-(4.45) and collecting the first order terms we have: s2 Jxh + hh s2I2 +12 s[Jih - J2/i] s2I2 + I2 = Ci >0 = c2 (4.104) (4.105) ICx + y/C'f + C2 2 - (4.106) (4.107) ~ //c,0 + s//c,l + S 2 / / c ,2 + - - -i i 2 ~ ao + sa\ + s a2... ~ ho + sh\ + s2h2 + 142 Chapter 4. Stability of multi layer flows 0(1) problem D(x'kDfkfi) - aoA,o = 0, A,o = 0 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,P{4>) + fkfi,h(4>), where Therefore, D(x'kDfkfi,P) = 0 /fc,o,p = 0 at 0 = 0,1, /fe,o,p = -*fe,o,0^o-/l,0,p = -/l0*2,0, 1 /2(0) '/2(0i)' (4.109) (4.110) where 1 Z"1 1 = / " 7 ^ J2(0) = / v Jo Xi J<t> X 2 d0. 143 Chapter 4. Stability of multi layer, hows Now, fko,h is a n eigen-function of D(x'kDfkfi,h) = <4fk,o,h, with OQ being the eigen-value. Therefore, from Sturm-Liouville theory, oft is real and less than zero. Thus, , ctR = 0 and aj 0. Our solution, fk,o{4>) = fkfi,p{4>) + /fe,o,/i(<^) has to satisfy the second jump condition, then ho f*W7S + *i,o,*7^) + X2Df2,o,h ~ x'iDflfi,h = -aih0t&nf3sm7r4>i, (4.111) equation (4.111) defines ho. Therefore, condition for instability is that A = — a2- is an eigen-value of both: D(x[Df)=\f /(O) = = o D(X2Df) = A/ fifo) = /(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 S t a t i c m u d c h a n n e l s As stated in Chapter 4, for certain 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 the flow, 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 | < 7 y / i ? . (5.3) tt((U,t) = o. (5.4) *(1 ,t,t) = Q(t) (5.5) *(0 ,O , t ) = * i „ f > , t ) , (5.6) *(4>,Z,t) =*mt(<t>,t), (5.7) (5.8) 145 Chapter 5. Static mud channels Whereas for the stability studies, we have Q(t) = 1, here we consider varying flow rate, such as would •be possible with pressure pulsing. For the velocity field we use an augmented Lagrangian algorithm, [26, 31], which fully describes the unyielded regions of the flow, the method is outlined in §5.1. We will work in a finite element approximation setting for the velocity field, with bilinear basis functions and square elements, for an summary of the finite element 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 the fluid is 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 a finite dimensional problem, i.e. since in any case we eventually discretise and solve a finite dimensional 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)=l-{Av,v)-{b,v) (5.9) 146 Chapter 5. Static mud channels where in (5.9), (•, •) is the canonical Euclidian inner product in R^. Let T be a linear mapping from RN into R M . Now, consider the minimisation problem J(u) < J(v) W e KerT (5.10) u e 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 {J(v) + (p,Tv)} (5-12) v e R 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 L(v,q) = J(v) + (q,Tv) (5.13) recall that (u, p) will be a saddle-point of L if L(u,q)<L(u,p)<L(v,p) Vv€R,N,qeRM (5.14) thus we have min max L(v, q) = max min L(v, q) = L(u, p) veRNqeR.M qeR.MveRN (5.15) It can be shown that L admits at least one saddle-point (u,p) on RN x R M , where u is the solution of (5.10) and is common to all the saddle-points of L on R w x R M . We define the augmented Lagrangian Lr , for r > 0, by Lr(v,q) = J(v) + (q,Tv) + ^\Tv\2 = L(v,q) + ^\Tv\2 (5.16) 147 Chapter 5. Static mud channels and note that a saddle-point of L(v,q) for which u G Ker T, will also be a saddle-point of Lr(v,q). However, (5.16) may be better conditioned to solve numerically. 5.1.2 A L G 2 In [26, 31] the authors describe an algorithm for finding the saddle-points of an augmented Lagrangian functional L r , they call it ALG2. Let us consider the minimization problem mm{F(Tv) + G(v)} (5.17) 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 = Vp and H = L2(il) x L2(ft) with their respective inner products and associated norms. We define, for v £ Vpfl, q € H and p G H the Lagrangian functional L(v,q,p)=F(q) + G(v) + (p,Tv-q) 148 Chapter 5. Static mud channels and then for r > 0, the augmented Lagrangian functional is defined as L{v,q,p)r = L(v,q,p) + r-\\Tv-q\\2. (5.19) The iterative method A L G 2 that we use below is in fact a method for calculating the saddle-points (u,p,X) of (5.19). A L G 2 . (p^A1) €H xH arbitrarily specified; (5.20) with (pn _ 1, A") known, determine succesively un,pn,Xn+1 by G(v)-G(un) + (Xn,T{v-un)) +r(Bun -pn~l,T(v - V)) > 0 Vv G Vp, (5.21) un G Vp F{q) - F{pn) - (An, q-pn) + r(pn - Tu11, q-pn))>0 (5.22) Va G H x H, and pn G H x H. yn+l =>n + pn(Tun _ pny (5.23) 5.1.3 Convergence of A L G 2 Results on convergence of A L G 2 are presented below, we just give a summary of what condi-tions 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) lim ^ = 0 0 (5.24) \Q\->°° \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: (Fo,(q)-F0,(p),q-p)>5M(\q-p\). (5.25) Now we are in a position to state the main result for convergence of the algorithm. Theorem 3 We assume that Lr admits a saddle-point (u,p,X) onVpxHxH. Under the assumption on T, F, G (i)-(iv) and if pn satisfies 0<Pn=P< (^Y^)r, (5-26) we have for ALG2 the following convergence results un —> u strongly in V, (5.27) pn —> p strongly in H, (5.28) 150 Chapter 5. Static mud channels Xn+l _ Xn _^ 0 strongiy i n H i (5 29) A" is bounded in H. (5.30) Proof For a proof see [26]. • 5.2 Application of A L G 2 to the displacement problem. 5.2.1 S t e a d y m o d e l Let us take T(u) = Vu, G(u) = 0, F(q) - F0(q) + Fi(q), where , /1 f\VV+q\* , 1/2) \ F0(q) = J UJ ^T-te + q-ndn, q € H x H, Fi(q) = I ^ | W * + g|dft, qeHxH. (5.31) 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 W1'1*1'™^) for $ and (L1'1+1/m(J7) x L^l+l/m{9)) for q. The 151 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 a finite-dimensional numerical method (Finite Element Method for example), the approximate solution space Vp^ is a finite dimensional 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 our finite dimensional solution space to be VPth C JY1(f2) and 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 fi • (Vv - q)dft + J / |Vv - q\2dfl. (5.32) Jn 1 Jn Therefore, if (u,p, A), where u e Vb,p and u is defined as * = + u, is a saddle point of L r , Vr > 0, then + 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 of finding un, pn, Xn+1 by solving sequentially: / A"-V(u-u") + (rVw"-p"-1)-V(7j-u")dft > 0, (5-33) Jn Vv G Vb,p, un e Vo,p F(q)-F(pn)+ f \ n • (q- pn) + r(pn — Vwn) • (q - pn)dfl > 0, (5.34) Jn Vq e H, pn e.H Xn+1 = \ n + P n ( V u n - p n ) , pn > 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, rAun = rV • pn~l - V • A n , (5 .36) where the right hand side is known from the previous step. The second step consists of mini-mizing equation (5 .34) with respect to q in order to obtain pn. This second step is the critical 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 .34) with respect to q in a discrete setting, we may locally minimize the integrand of Lr(un, q, A " ) with respect to q G Hh x Hh where Hh denotes the restriction of if to a particular discrete element, thus pn = mi{-Jo ^ L 7 - i d s + g - f + ^ | V * * + q | + - | q | 2 - ( A " + r V U " ) - q } , (5 .37) equivalently we have, 1 /.|v**+g| 2 / 1/2^ p n = lqi{2 j ^^ds+^ |V**+g |+ - |V**+q | 2 - (A"+rV(r+«") - / ) - (V**+ g ) } . (5 .38) Clearly, this expression will achieve its minimum when (V** + q) || ( A n + rV(** + un) - f). Therefore, letting ( V * * + p " ) = 0 ( A " + rV(** + u")-/), (5 .39) x = |A n +rV(**+u")-/ | (5 .40) we have to find the minimizer of 1 r(0x)2 Y ( S J / 2 ) T V r o M { 0 ) = 2jo ^ 7 ^ d s + + 2 ( 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, X(9x) + ^ F + r 0 x - x = 0, • ( 5 . 4 2 ) ' / H which has a unique solution and therefore here the fluid is yielded in this iteration. Finally we update A using ( 5 . 3 5 ) . 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, Lr(vn+\qn+\n) = F(qn+1) + f pi • ( W * + 1 - q n + 1 ) d f 2 + J / |Vi> n + 1 - q" + 1 | 2 df2 Jn * Jn (5.43) where the super index vn+1 means v(tn+l) where tn+1 =tn + At. Thus the transient version of A L G 2 will be, jT \ n k + l • V(vn+1 - unk+l) + ((r +• £)Vu£+1 - -LVttf - pnkt\) • V(«" + 1 - x)dQ > ' 0, Vvn+1 G Vo,p, unk+1 e V0,p (5.44) F(qn+1) - F(p£ + 1 ) + f \n+l • (qn+1 - pnk+1) + r(p£ + 1 - qn+1) • (qn+1 - pJJ+ 1)dfi > 0, Jn V q " + 1 G H, pl+1 G H (5:45) KX\=K+1+P^K+1-Pl+l), 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 < + 1 = r V - p n - 1 - V - A " - ^ V ^ . (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,p such that a(u,v) = b(y) V v £ Vo,p (5.48) where a(u, v) = / Vu • Vudfi, • Jn b(v) = [ fvdtl, Jn with / being the right hand side of (5.36). The finite element discretization consists in replacing VQIP by a finite dimensional subspace Vh, i.e. instead of solving (5.48), we solve Find uh £ Vh such that a(uh,v) = b(v) Vv£Vh, (5.49) 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' £Th, K ^ K' and interior^) f| interior '^) = 0. 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 Th}, 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, Vh = {u £ C(U)\ u\K £ P\{K) \/K £ Th, u = 0 on dfl,} . 155 Chapter 5. Static mud channels where P\ (K) is the set of polynomials of first degree (in 2 variables), i.e. pG P\{K) & p(x,y) = a + Bx + 7y 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, a2, 0 3 are the vertices of K, thus p(x, y) = ui<pi(x, y) + u2if2(x, y) + U3ip3(x, 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 Tn. 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 a finite volume 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 disper-sion, 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 the flux by 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 numerical fluxes across 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 transport fluxes across each element and using a forward Euler scheme 157 Chapter 5. Static mud channels for time we have: /~m+l /-m At , „ n j - , n s At, n n . (5.52) where F and G are the discrete numerical fluxes. Note that (5.52) is a second order approx-imation 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 FL r . and GL 1 using a low order scheme in the following way, and G Vni .CJ 1 , HVN, . >0 I V7! x .CS.!,- i f V " ! . <0 w . j . i C j 1 , - i f i y p . ^ j ^ o r . . c ^ i if w ™ , < o . Now we compute the low order time advanced solution, A*-£ F 1 , .) — x 7 ( C f . , i — GL. 0-?.— ± . 1 ' \ 7.1-1 -4 *).? — ~ AC Then we compute F^f i . and G f f. i using a high order scheme in the following way (5.53) (5.54) (5.55) (5.56) (5.57) We compute the antidiffusive fluxes, to minimize the diffusion introduced by using a low order scheme, tun = K+h+Gi^-(5.58) (5.59) We limit this antidiffusive fluxes to minimize the dispersion effect that the high order approxi-mation has on the solution, A1. (5.60) (5.61) 158 Chapter 5. Static mud channels and finally we apply the limited antidiffusive fluxes to the get the final solution, - ( (5.62) The only thing that remains unclear is how the corrections Li+i j and Lt.j+i will be chosen to limit C^f1 so that it does not exceed a maxim value Cj"?x nor a minimum value CJfn. This will be achieved in the following way. First we define the following quantities, 1. We compute the sum of all antidiffusive fluxes into the grid point (i, j) (the centre of our element), P± = max(0,^_i j) - mm(0,Ai+ij) + max(0, ) - 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, Aid_i). 3. Define 4. Define s-rmax _ m a v ( r < n rm rm rm rm \ i,j — m m V ° i , j ' ° i - l j > (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 (5.63) (5.64) 159 Chapter 5. Static mud channels 6. Finally, the corrections become min(i?+.+ 1,i? i j) min(JR+.,i?- j+1) nan(Rt+1j,Rid) min(i?+.,i?-+lj.) if A-, if A - , i f V * J < 0 (5.65) (5.66) Now, we are in the position to solve the system of equations (5.1)-(5.7) numerically and inves-tigate 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 pulsating flow rate, thus (5.5) becomes Thus we will solve system (5.1)-(5.7) with (5.67) instead of (5.5), using the augmented La-grangian algorithm in a finite element setting for the velocity field and 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 to find an analytical solution that de-scribes 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 and final similar 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. 5P = 0 in (5.67). The following physical and rheological parameters are used, K\ = 0.5, K2 = 0.4, m\ = 1, m2 = 1.2, px = i ) p2 = 0.9, 7V,i = 1, 7V,i = 0.9, 7Y , ! = 0.7 e = 0.0 and 3 = 0. Results are consistent #(!,£, t) = 1 + dpsincot, (5.67) 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 (5P = 0). 0.8 0.6 0.4 0.2 N 0 -0.2 -0.4 -0.6 -0.8 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, K2 = 0.4, m\ = 1, m2 = 1.2, P l = l ; P 2 = 0.9, 7y , i = 1, ry,i = 0.9, TY,I = 0.7 e = 0.0 and (3 = 0. 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, Ty,i = 1, ryti = 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 0 t t (a) (b) (c) 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 = 2-K/U), for the transient velocity model is defined as T = 2TT./U)€. 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 5P = 0.1. No major effects can be seen on the displacement front and 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 the flow rate 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, Sp = 0.1. As 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) both fluids are 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 displacing fluid is forced to the narrow side and the displaced fluid is 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 8P = 0.1 and 5P = 0.2 respectively. As in the previous cases there is no clear evidence of a change in the displacement. Figures 5.16 and 5.17 show the effects of a pulsating flow rate on the transient model with amplitudes, Sp = 0.1 and 5P = 0.2 respectively. Clearly the effects of 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 (a) (b) (c) (d) 4> <f $ ' • (e) (f) (g) (h) Figure 5.3: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Steady state velocity model. Period T = 2TX/U), U = 10 and Sp = 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, m\ = 1, m2 = 1.2, p\ = 1, p2 = 0.9, Ty,i = 1, 7 y , i = 0.9, Ty,i = 0.7 e = 0.3, /3 = 0. 165 Chapter 5. Static mud channels (a) (b) (c) (d) in HI flMI 0.2 0.5 i i ' " 0.8 (e) (f) (h) Figure 5.4: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Steady state velocity model. Period T = 2n/uj, OJ = 10 and Sp = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 — 0.9, Ty,\ — 1, ry,! = 0.9, 7 y i = 0.7 e = 0.3, /3 = 0. 166 Chapter 5. Static mud channels (a) (b) (c) (d) (e) (f) (g) (h) Figure 5.5: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Transient velocity model. Period T = 27T/W, UJ = 10, 5V = 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: KI = 0.5, « 2 = 0.4, mi = 1, m2 = 1.2, p\ = I, p2 — 0.9, Ty,i = 1, TY,I = 0.9, 7y , i = 0.7 e = 0.3, /3 = 0. 167 Chapter 5. Static mud channels (a) (b) (c) (d) (e) . (f) (g) (h) Figure 5.6: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Transient velocity model. Period T = 2-n/uie, u> = 10, 5P = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, ry^ = 1, ry,i = 0.9, Ty,i = 0.7 e = 0.3, 8 = 0. 168 Chapter 5. Static mud channels 0.2 0.5 0.8 0.2 0.5 (a) (b) (c) (d) (e) (f) (g) (h) Figure 5.7: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Transient velocity model. Period T = 2n/uje, UJ = 10, Sp = 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, K2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, Tyt\ = 1, TY,I = 0.9, TY,I = 0.7 e = 0.3, /3 = 0. 169 Chapter 5. Static mud channels m 3.5 3 2.5 2 1.5 1 0.5 0.2 0.5 0.2 0.5 0.8 (a) (b) (c) (d) 4HH r 2.5! - \ 0.2 0.5 (e) (f) (g ) (h) Figure 5.8: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Transient velocity model. Period T = 2n/u>e, OJ = 1, Sp = 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, m2 = 1.2, p\ = 1, p2 = 0.9, Ty,i = 1, TYI = 0.9, TYi = 0.7 e = 0.3, 8 = 0. 170 Chapter 5. Static mud channels 0.2 0.5 0.8 0.2 0.5 0.8 4> <f . (a) (b) Figure 5.9: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Transient velocity model. Period T = 27r/we, UJ = 0.1, 5p = 0.1 and e = 10. Interface position and velocity field for times: a)-b) T/4. Physical and rheological parameters: K\ = 0.5, K 2 = 0.4, mi = 1, m2 = 1.2, P l = 1, p2 = 0.9, TyA = 1, TY,I = 0.9, 7y , i = 0.7 e = 0.3, /5 = 0. 171 Chapter 5. Static mud channels (a) (b) (c) (d) <)> <> «t> <t> - • (e) (f) (g) (h) Figure 5.10: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Steady state velocity model. Period T = 2-ir/uj, UJ = 10 and 5P = 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 2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, Ty,i = 1, TY,I = 0.9, ryA = 0.7 e = 0.5, (3 = 0. 172 Chapter 5. Static mud channels J m : III iii!: •ill: 0.5 0.8 • ' I n In In In ii iiii n, I . , . l i t 0.2 0.5 0.8 (a) (b) (c) (d) ! ! : i n i n " i n in •i! I l l III III 0.2 0.5 0.8 i i . - i i i in In III 1 III,, i i , ; ; (e) (f) (h) Figure 5.11: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Steady state velocity model. Period T = 2TT/UJ, U> = 10 and 5p = 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, K2 = 0.4, mi = 1, m2 — 1.2, p\ = I, p2 = 0.9, ry,\ = 1, Ty , i = 0.9, 7 y i = 0.7 e = 0.5, 8 = 0. 173 Chapter 5. Static mud channels o o o o (a) (b) (c) (d) 0 <f «. 0 « (e) (f) (g) (h) Figure 5.12: Displacement flow in eccentric annulus, interface propagation. Unsteady displace-ment. Transient velocity model. Period T = 2ir/ue, OJ = 10, 5p = 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: KI = 0.5, «2 = 0.4, mi = 1, TO2 = 1.2, p\ = 1, p2 = 0.9, ry,i = 1, TY,I = 0.9, ry,! = 0.7 e = 0.5, (3 = 0. 174 Chapter 5. Static mud channels Figure 5.13: Displacement flow in eccentric annulus, interface propagation. Mud channel for-mation. Steady state velocity model. Period T = 2TT/U, to = 10 and Sp = 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, m2 = 1.2, p\ = 1, p2 = 0.9, Ty,i = 1, ry,i = 0.9, TY,I = 0.7 e = 0.8, /3 = 0. 175 Chapter 5. Static mud channels 41 0.2 0.5 ,-»•. J 0.2 0.5 0.8 (a) (b) (c) (d) ""•Mil i, i i,i in in 1 In in in i n 05 • n iii , , , , 0.2 0.5 0.8 0.2 0.5 (f) (h) Figure 5.14: Displacement flow in eccentric annulus, interface propagation. Mud channel for-mation. Steady state velocity model. Period T = 2n/u, u = 10 and 6P = 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 2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p 2 = 0.9, Ty,i = 1, Ty,i = 0.9, 7 v , i = 0.7 e = 0.8, /3 = 0. 176 Chapter 5. Static mud channels ^ 1.5 2.5fc *M.5 0.5 [ 0.2 0.5 0.8 ^ 1.5 0.2 0.5 0.8 i n ' 2.5:; . . . 0.5h *1 0.2 0.5 0.8 (a) (b) (c) (d) 1.5 i n * 2.5;, I •1.5 0.5 0' In 0.2 0.5 0.8 "J-1.5 2.5 In In j ! 0.5 0.2 0.5 0.8 (e) (f) (g) (h) Figure 5.15: Displacement flow in eccentric annulus, interface propagation. Mud channel for-mation. Steady state velocity model. Period T = 2-rr/uj, LO = 10 and 5p = 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. 177 Chapter 5. Static mud channels 8811 0.2 0.5 0.8 0.2 0.5 0.8 0.2 0.5 0.8 (a) (b) (c) (d) 0.2 0.5 0.8 0.2 0.5 0.8 (e) (f) (g) (h) Figure 5.16: Displacement flow in eccentric annulus, interface propagation. Mud channel for-mation. Transient velocity model. Period T = 2ir/coe, u> = 10, 5P = 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 2 = 0.4, mi = 1, m2 = 1.2, pi = 1, p2 = 0.9, TY,I = 1, 7y, i = 0.9, 7y, i =0.7 e = 0.8, 8 = 0. 178 Chapter 5. Static mud channels (a) (b) (c) (d) • <t> • • (e) (f) (g) (h) Figure 5.17: Displacement flow in eccentric annulus, interface propagation. Mud channel for-mation. Transient velocity model. Period T = 2-jr/coe, co = 10, Sp = 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. 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 Sp = 0.2. This model predicts that 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 5P = 0.2. Clearly as we go over the full period of 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 Sil l l l l l l i i l ikE •••((••In • H I • S H ^ ^ ^ ^ ^ 0.2 0.5 0.8 (a) (b) (c) (d) (e) (0 Figure 5.18: Displacement flow in eccentric annulus, interface propagation. Mud channel. Steady state velocity model. Period T = 2-ir/cj, ui = 10 and 5P = 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 2 = 0.4, mi = 1, m2 = 1.2, pi = 1, p2 = 0.9, Ty,i — 1, r^i = 0.9, = 0.7 e = 0.8, 8 = 0. 181 Chapter 5. Static mud channels (a) (b) (c) (d) (e) ' (f) (g) (h) Figure 5.19: Displacement flow in eccentric annulus, interface propagation. Mud channel. Transient velocity model. Period T = 2jr/ue, u = 10, 5P = 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, n2 = 0.4, mi = 1, m2 = 1.2, p\ = 1, p2 = 0.9, ry,\ = 1, 7y,i = 0.9, Tyti = 0.7 e = 0.8, 3 = 0. 182 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 demon-strated 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 long finger develops we find essentially a multi-layer parallel flow along the length of the annulus. In chapter 4 we studied the linear stability of these flows for the case when both fluids are 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, parallel flow of two Newtonian fluids and parallel flow of 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 the flow is linearly stable. Inclination of the Hele-Shaw cell has destabilizing effects, if the heavier fluid is on top of the lighter, fluid, the flow may 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 stress fluids has to be solved via computation. Here we solve the full problem, i.e. the parallel flow of two yield stress fluids in a vertical or near vertical eccentric Hele-Shaw cell. We consider the case when both fluids are yielded at the interface. If the eccentricity is equal to zero, i.e. a concentric annulus, we find analogous results to the Newtonian and power law fluid cases: 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 the flow may be stable, but if p\ > pi and n\ < K2 the flow may be 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 2 , the flow was linearly unstable for all wave lengths. If we increase 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 the flow will 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 the flow will 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 parallel flow may be unstable. Therefore, in unsteady displace-ments 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 parame-ter regimes for which the annular displacement would be steady. Note that for such parameters, a base parallel flow is a perfectly well-defined solution of the two fluid problem. 6.3 Mud channel removal In chapter 5 we studied the effects of a pulsatile flow rate for the removal of a mud channel. This is via numerical simulation. In chapter 4 we showed that if one of the fluids is unyielded at the interface the flow will be linearly stable. We investigated the effects that pulsation has on the mud channel by first using 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 the flow rate. 186 B i b l i o g r a p h y F . ANDREU, C. BALLESTER, V. CASELLES, J.M. MAZON, The Dirichlet Problem for the Total Variation Flow, Journal Functional Analysis,180, pp. 347-403, (2001). R. ARIS, On the Dispersion of a Solute in a Fluid Flowing through a Tube, Proc. Roy. Soc. A (235), pp. 67-77, (1956). R. ARIS, On the Dispersion of a Solute in Pulsating Flow Trhough a Tube, Proc. Roy. Soc. A (259), pp. 370-376, (1960). J. AZAIEZ A N D B. SINGH, Stability of miscible displacements of shear thinning fluids in a Hele-Shaw cell. Phys. Fluids. (14)5, pp. 1557-1571, (2002). R. BALASUBRAMANIAM, N . RASHIDNIA, T. MAXWORTHY A N D J. KUANG, Instability of miscible interfaces in a cylindrical tube . Phys. Fluids 17, 052103, pp. 1-11, (2005). G. I. BARENBLATT, V.M. ENTOV A N D V.M. RYZHIK, Theory of Fluid Flows Through Natural Rocks. Kluwer Academic Publishers (1990). H . A. BARNES, J.F. HUTTON A N D K. WALTER, An introduction to Rheology. Elsevier (1989). G.K. BATCHELOR, An Introduction to Fluid Dynamics. Cambridge University Press (1989). R.M. BEIRUTE A N D R.W. FLUMERFELT, Mechanics of the Displacement Process of Drilling Muds by Cement Slurries Using an Accurate Rheological Model, 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 A N D I.A. F R I G A A R D , Mud removal and cement place-ment during primary cementing of an oil well; laminar non-Newtonian displacements in an eccentric annular Hele-Shaw cell. J. Engng. Math., 43, pp. 229-253, (2002). J.W. B R I C E A N D R.C. H O L M E S , Engineered Casing Cementing Programs Using Turbulent Flow Techniques, J. Pet. Tech., pp. 503-508 (1964). M.A. CHILDERS, Primary Cementing of Multiple Casing, Society of Petroleum Engineers, paper number SPE 1919, (1968). C. R. CLARK A N D L . G . CARTER, Mud Displacement with Cement Slurries, Society of Petroleum Engineers, paper number SPE 4090, (1973). M . COUTURIER, D . GUILLOT, H . HENDRIKS A N D F. CALLET, Design Rules and Associ-ated Spacer Properties for Optimal Mud Removal in Eccentric Annuli, Society of Petroleum Engineers, paper number SPE 21594 (1990). P. COUSSOT, Mudflow Rheology and Dynamics, A.A. Balkema Rotterdam (1997). P. COUSSOT, Saffman-Taylor instability in yield-stress fluids. J. Fluid Mech. (380), pp. 363-376, (1999). 187 Bibliography [17] B. D A C O R O G N A , Direct Methods in the. Calculus of Variations. AMS Vol. Springer-Verlag (1989). [18] P.G. DRAZIN, Introduction to Hydrodynamic Stability. Cambridge Univeristy Press (2002). [19] D . DUSTREHOFT, G . WILSON A N D K . NEWMAN, Field Study on the Use of Cement Pulsation to Control Gas Migration. Society of Petroleum Engineers paper number SPE 75689 (2002). [20] G . DUVAT A N D J.L. LIONS, Inequalities in mechanics and physics. Springer-Verlag (1976). [21] M.J. ECONOMIDES, IMPLICATIONS O P CEMENTING O N W E L L PERFORMANCE. IN: 'E. B. NELSON ( E D . ) , Well Cementing, Schlumberger Educational Services (1990) Chapter 5. [22] I. EKELAND, R. TEMAM, Convex Analysis and Variational Problems. SIAM (1999). [23] P. ERN, F. CHARRU A N D P. LUCHINI, Stability analysis of a shear flow with strongly stratified viscosity. J. Fluid Mech. 496, pp.295-312, (2003). [24] L.C. EVANS, Partial Differential Equations. American Mathematical Society, (1998). [25] J. FERNANDEZ, P. KUROWSKI, P. PETITJEANS A N D E. MEIBURG, Desntiy-driven unsta-ble flows of miscible fluids in a Hele-Shaw cell. J. Fluid Mech. 451, pp. 239-260, (2002). [26] M. FORTIN A N D R. GLOWINSKI, Augmented Lagrangian Methods. North-Holland, (1983). [27] LA. F R I G A A R D , S.D. HowiSON A N D I.J. SOBEY, On the Stability of Poiseuille Flow of a Bingham Fluid. J. Fluid Mech. 263, pp. 133-150, (1994). [28] I.A. FRIGAARD, Super-stable parallel flows of multiple viscoplastic fluid. J. Non-Newt. Fluid Mech., 100, pp 49-76, (2001). [29] D . GILBARG, N. S. TRUDINGER, Elliptic partial differential equations of second order. Berlin ; New York : Springer-Verlag, (1983). [30] R. GLOWINSKI, J.L. LIONS A N D R. TREMOLIERES, Numerical analysis of variational inequalities. North-Holland (1981). [31] R. GLOWINSKI, Augmented Lagrangian Methods: Applications to the numerical solution of boundary-value problems. North-Holland (1983). [32] P. GONDRET, N. RAKOTOMALALA, M. RABAUD, D . SALIN A N D P. WATZKY, Viscous parallel flows in finite aspect ratio Hele-Shaw cell: Analytical and numerical results. Phys. Fluids 9(6) (1997). [33] P. GONDRET A N D M. RABAUD, Shear,instabillity of two-fluid parallel flow in a Helle-Shaw cell. Phys. Fluids 9(11) (1997). [34] V.A. GORODTSOV A N D V.M. YENTOV, Instability of the displacement fronts of non-Newtoinan fluids in a Hele-Shaw cell. J. Appl. Maths. Mechs. (61)1, pp. 111-126 (1997). . [35] R. GOVINDARAJAN, Effect of miscibility on the linear instability of two-fluid channel flow. Int. J. Multiphase Flow. (30), pp. 1177-1192, (2004). 188 Bibliography [36] N. G O Y A L A N D E. M E I B U R G , Unstable density stratification of miscible fluids in a vertical Hele-Shaw cell: influence of variable viscosity on the linear stability. J . Fluid Mech. 516, pp. 211-238, (2004). [37] J . H A R R I S , Rheology and Non-Newtonian flow. Longman (1977). [38] E . J . H I N C H , A note on the mechanism of the instability at the interface between two shearing fluids. J . Fluid Mech., 114, pp. 463-465, (1984). [39] E . J . H I N C H A N D F. P L O U R A B O U E , Kelvin-Helmholtz instability in a Hele-Shaw cell: Large effect from the small region near the meniscus. Phys. Fluids. (17)052107, pp. 1-13, (2005). [40] G . M . H O M S Y , Viscous Fingering in porous media. Ann. Rev. Fluid Mech. (19), pp. 271-311, (1987). [41] A.P. H O O P E R A N D W . G . B O Y D , Shear flow instability at the interface between two viscous fluids. J . Fluid Mech., 179, pp. 201-225, (1987). [42] 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 Methods for Elliptic and Parabolic Partial Differential Equations. Springer-Verlag New York (2003). [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 of non-Newtonian Hele-Shaw flow. Physical Rev. E. (54)5, pp. R4536-R4539, (1996). [45] 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 Hele-Shaw flow and the Saffaman-Taylor instability. Physical Rev. Lett. (80)7, pp. 1433-1436, (1998). [46] 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. S A L I N , 3D Instability of Mis-cilbe Displacements in a Hele-Shaw Cell. Phys. Rev. Lett. (79)26, pp. 5254-5257, (1997). [47] 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 , D. S A L I N A N D Y.C. Y O R T S O S , Mis-cible displacement in a Hele-Shaw eel at high rates. J . 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 , Viscous fingering in a yield stress fluid. Phys. Rev. Lett. (85)2, pp. 314-317, (2000). [49] H.P. L A N G T A N G E N , Computational Partial Differential Equations Springer-Verlag Berlin Heidelberg (1999). [50] C F . L O C K Y E A R A N D A.P. H I B B E R T , Integrated Primary Cementing Study Defines Key Factors for Field Succes, J . Pet. Tech. pp. 1320-1325, (1984). [51] C F . L O C K Y E A R , D.F. R Y A N A N D M.M. G U N N I N G H A M , Cement Channeling: How to Predict and Prevent, Society of Petroleum Engineers, paper number SPE 19865, (1990). [52] R.H. 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 Mechanics in Pri-mary Cementing, J . Pet. Tech. pp. 251-260, (1967). 189 Bibliography [53] H . P A S C A L , Rheological behaviour effect of non-Newtonian fluids on dynamics of moving interface in porous media. Int. J. Engng. Sci. (22), pp. 227-241 (1984). [54] H . P A S C A L , Dynamics of moving interface in porous media for power law fluids with a yield stress. Int. J. Engng. Sci. (22), pp. 577-590 (1984). [55] H . P A S C A L , A theoretical analysis of stability of a moving interface in a porous medium for Bingham displacing fluids and its application in oil displacement 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 instability of long fingers in Hele-Shaw flows Phys. Fluids. (28)6, pp. 1583-1585 (1985). [57] 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 in primary cementing of an oil well. J. Engng. Math., 46(1), pp.1-26, (2004). [58] S. P E L I P E N K O & I. A. F R I G A A R D , Two-dimensional computational simulation of eccentric annular cementing displacements. IMA J. Appl. Math., 69, pp. 557-583, 2004. [59] S. P E L I P E N K O & I.A. F R I G A A R D , Visco-plastic fluid displacements in near-vertical nar-row eccentric annuli: prediction of travelling wave solutions and interfacial 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 , Kelvin-Helmholtz instability in a Hele-Shaw cell. Phys. Fluids. (14)3, pp. 922-929, (2002). [61] R. R A G H A V A N A N D S.S. M A R S D E N , The stability of immiscible liquid layers in a porous . medium. J. Fluid Mech. (48), pp. 143-159, (1971). [62] R. R A G H A V A N A N D S.S. M A R S D E N , A theoretical study of the instability in the ^parallel flow of immiscible liquids in a porous medium. 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 Techinque, Society 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 , The penetration of a fluid into a porous medium or Hele-Shaw cell containing a mor viscous liquid, 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 , Stability and Transition in Shear Flows. Springer-Verlag, New York, Inc., (2001). [67] 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 instabilities during displace-ments of two miscible fluids in a vertical pipe . Phys. Fluids. (13)3, pp. 553-556, (2001). [68] M. S H A R I A T I A N D Y.C. Y O R T S O S , Stability of miscible displacements across stratified porous media. Phys. Fluids. (13)8, pp. 2245-2257 (2001). 190 Bibliography [69] R.C. S I M I T H , Succesful Primary Cementing Can Be a Reality, J . Pet. Tech. pp. 1851-1858, (1984). [70] T.R. S M I T H , Cementing Displacement Practices-Field Applications, J . Pet. Tech. pp. 564-629, (1984). [71] J.P. S T O K E S , D . A . W E I T Z , J.P. G O L L U B , A . D O U G H E R T Y , M.O. R O B B I N S , P.M. C H A I K I N A N D H . M . L I N D S A Y , Interfacial stability of immiscible displacement in a porous medium Physical Rev. Lett. (57)14, pp. 1718-1721 (1986). [72] P. S Z A B O A N D O. H A S S A G E R , Flow of Viscoplastic Fluids in Eccentric Annular Geome-tries. J . Non-Newt. Fluid Mech, 45, pp 149-169, (1992). [73] P. S Z A B O A N D O. H A S S A G E R , Displacement of One Newtonian Fluid by Another: Denstiy Effects in Axial Annular Flow, Int. J . Multiphase Flow. (23)1, pp. 113-129, (1996). [74] C.T. T A N A N D G . M . H O M S Y , Stability of miscible displacements in porous media: Recti-linear flow. Phys. Fluids. (29)11, pp. 3549-3556 (1986). [75] C.T. T A N A N D G . M . H O M S Y , Simulation of nonlinear viscous fingering in miscible dis-placement. Phys. Fluids. (31)6, pp. 1330-1338 (1988). [76] G . T A Y L O R , Dispersion of Soluble Matter in Solvent Flowing Slowly through a Tube, Proc. Roy. Soc. A (219), pp. 186-203, (1953). [77] G . T A Y L O R , Conditions under Which Dispersion of a Solute in a Stream of Solvent can be Used to Measure Molecular Diffusion, Proc. Roy. Soc. A (225), pp. 473-477, (1954). [78] A . T E H R A N I , S .H. B I T T L E S T O N & P.J. L O N G , Flow instabilities during annular dis-placement of one non-Newtonian fluid by another. 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 stability theory of two-layer fluid flow in a inclined channel. Phys. Fluids. (6)12, pp. 3906-3922, (1994). [80] 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 unsteady flows of Bingham fluis: a numerical strategy and some benchmark results. J . Comp. Phys. (187), pp. 441-456, (2003). [81] D . V O L A , F. B A B I K A N D J . C . L A T C H E , On a numerical strategy to compute gravity currents of non-Newtonian fluids. J . Comp. Phys. (201), pp. 397-420, (2004). [82] G . V I N A Y , A . W A C H S A N D J .F . A G A S S A N T , Numercial simulation of non-isothermal viscoplastic waxy crude oil flows. J . Non-Newt. Fluid Mech, 128, pp. 144-162, (2005). [83] 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 plastic in a narrow eccentric annulus. J . Fluid Mech. 222, pp. 39-60, (1991). [84] S.D.R. W I L S O N , The Taylor-Saffman problem for a non-Newtonian liquid. J . Fluid Mech. (220), pp. 413 425, (1990). [85] Z. Y A N G A N D Y . C . Y O R T S O S , Asymptotic solutions of miscible displacements in geome-tries of large aspect ratio. Phys. Fluids. (9)2, pp.286-298, (1997). 191 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 regimes in unstable miscible dis-placements in random porous media. Advances in Water Resources. (25), pp. 885-898, (2002). [87] S.G. Y I A N T S O S A N D B.G. HlGGINS, Linear stability of plane Poiseuille flow of two su-perposed fluids. Phys. Fluids, (31) ,pp. 3225-3238 (1988). [88] C.-S. Y I H , Instability due to viscosity stratification. J. Fluid Mech, 27(2), pp. 337-352, (1967). [89] S.T. Z A L E S A K , Fully Multidimensional Flux-Corrected Transport Algorithms for Fluids. J. Non-Newt. Fluid Mech, (100), pp 49-76, (2001). [90] E. Z E I D L E R , Applied functional analysis : applications to mathematical physics. Springer-Verlag, (1995). [91] M. Z E Y B E K A N D Y . C . Y O R T S O S , Parallel flow in Hele-Shaw cells. J. Fluid Mech. (244), pp. 421-442, (1992). [92] J. Z H A N G A N D LA. F R I G A A R D , Dispersion effects in the miscible displacement of two fluids in a duct of large aspect ratio. 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 Cementation, Society of Petroleum .Engineers, paper number SPE 4830, (1975). 192
- 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
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
Page Metadata
Item Metadata
Title | Transient effects in oilfield cementing flows |
Creator |
Moyers González, Miguel Angel |
Publisher | University of British Columbia |
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 |
GraduationDate | 2006-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | 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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0080123/manifest