NONLINEARLY STABLE MULTILAYER VISCOPLASTIC FLOWS. by MIGUEL A N G E L MOYERS G O N Z A L E Z BSc (Applied Mathematics) I .T.A.M., 2000 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA August 2002 © Miguel Angel Moyers Gonzalez, 2002 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Date /Q - Sep f - ZOOZ Abstract When two purely viscous fluids flow in a parallel multi-layer shear flow down a pipe or other duct, the flow is frequently unstable. The primary instability arises at the interface, which is linearly unstable to certain wavelengths. Physically, for any perturbation of a planar interface, the less viscous fluid accelerates the more viscous fluid; see [5]. Although it is possible to suppress short wavelength instabilities through surface tension and also to achieve a more stable flow through arranging the viscosities of the fluid layers, (see e.g. [9]), none of these methods actually eliminate the underlying interfacial instability. In this thesis we demonstrate that by using a visco-plastic fluid, (i.e. a fluid with a yield stress), as the lubricating fluid in a multi-layer pipe flow, we do actually eliminate interfacial instabilities. This fact has been established recently in [2] using a linear stability analysis. In this thesis the analysis is extended to nonlinear stability, applying the methods developed in [10]. Our basic flow is an axisymmetric flow in which a Bingham fluid surrounds a Newtonian fluid. The flow rates and rheologies are such that the Bingham fluid is unyielded at the interface. We consider the nonlinear stability of this flow. This nonlinear stability problem is interesting for a number of reasons. Firstly, for a nonlinear perturbation it is possible for the unyielded plug region, (and the New-tonian region that it surrounds), to be perturbed and move. For the linear stability problem in [2], such motions do not occur since the stress perturbations on the interface and yield surface are linear and periodic, (via a normal mode approach). Thus, they integrate to zero over the volume of the plug and there is no motion. To allow for motion of the plug and Newtonian fluid region, whilst preserving a stable interface, is the main challenge of this thesis. The stable interface is preserved by bounding the perturbation in the deviatoric stress, which ensures that an unyielded region will still surrounds the Newtonian region. To carry out the analysis for a moving region we then consider a perturbation about a basic flow that is asym-metric. We find stability bounds, (conditional on the stress perturbation), below which the velocity perturbation from the asymmetric basic flow decays. We then show that the displace-ment of the plug region is also bounded for all time. Thus, in the above way we prove that the axisymmetric basic flow is stable but not asymptotically stable, i.e. conditional nonlinear perturbations will decay to an asymmetric solution that is close to the axisymmetric basic flow. Secondly, the asymmetric basic flows are themselves interesting. These flows fall into a class of flows that have been studied in [3], from the point of view of proving existence and uniqueness, in the context of examining various problems in cementing. Here we compute these flows using the augmented Lagrangian approach. This method is special in that it manages to compute the unyielded behaviour exactly in regions where the stress is below the yield stress. This is the first time, to our knowledge, that this method has been used to compute multi-fluid flows of visco-plastic fluids. Thirdly, the flows we study are special in that the structure is preserved, (i.e. due to one of the fluids having a yield stress). For any multi-phase fluid problem this is unusual. It is really this preservation of structure that has allowed the application of nonlinear stability methods. As such, this study is one of very few multi-phase flow problems that we know of, where an energy method has been effectively applied. ii Abstract Finally, in demonstrating weak nonlinear stability of these flows, we open the way for a range of applications to industrial processes where multi-layer processing might be of interest, e.g. coated wires and tapes, co-extruded laminates, etc. iii Table of Contents Abstract " Table of Contents iv List of Tables vi List of Figures vii Acknowledgements viii Chapter 1. Introduction 1 1.1 Problem description 1 1.1.1 Applications 1 1.2 Background literature 2 1.3 Outline of the thesis 4 Chapter 2. Viscoplastically lubricated flows 6 2.1 Model derivation 7 2.1.1 Equations of motion 7 2.1.2 Dimensionless Equations of motion 8 2.2 Axisymmetric Basic flows 10 2.2.1 Determingr* 14 Chapter 3. Numerical solution of the asymmetric basic flow 18 3.1 Asymmetric basic flows 19 3.1.1 Variational form of the problem 20 3.2 Augmented Lagrangian Method 22 3.2.1 Principles of the method 23 3.2.2 Algorithms using Augmented Lagrangian Methods 24 3.2.3 Application of A L G 2 to our problem 33 3.3 Example Results 35 3.3.1 Convergence 36 3.3.2 Preservation of the plug region 37 3.3.3 Different solution types for asymmetric basic flows 41 3.3.4 Variations with B 41 3.3.5 Variations with m 44 3.3.6 Variations with Ri 44 Chapter 4. The Reynolds-Orr equation and Stability bounds 48 4.1 Derivation of the Reynolds-Orr equation 49 4.2 Three-dimensional nonlinear stability bounds 51 4.2.1 The Bingham fluid region, fi2(i) 5 2 4.2.2 The Newtonian fluid region, fii(i) 56 iv Table of Contents 4.3 The stability bound 61 4.3.1 The bound for \{xc,yc)\ 62 4.4 Parametric variations in the stability bounds 63 4.4.1 \VW\2,max, \VcW\max, WN,max 63 4.4.2 AB,i, A B ) 2 64 4.4.3 AJV,I, Aiv,2, AAT,3, AAT,4 65 4.4.4 Variation of (4.49) & (4.50) 67 Chapter 5. Discussion and conclusions 68 5.1 The bounds (4.49) & (4.50) 68 5.1.1 Stability where |(:r c,y c)| = 0 69 5.1.2 One-dimensional stability 70 5.2 Summary of the main contributions of the thesis 70 5.3 Future work 70 Bibliography 72 v List of Tables 3.1 Comparison of various asymmetric basic flows 3.2 Variations with Bingham number. 3.3 Variations with viscosity ratio, m. 3.4 Variations with i?, List of Figures 2.1 (a) Newtonian fluid surrounded by a Bingham fluid flowing in a pipe , (b) cross section of the pipe fii corresponds to the Newtonian fluid, fi2 corresponds to the Bingham fluid, (c) Desired velocity profile (Newtonian fluid surrounded by an unyielded plug region of Bingham fluid) 7 2.2 Basic flow for Ri < r* < 1 (Ri = 0.5, m = 38 and B = 100) 13 2.3 Basic flow for r* < Ri and B < p 4 , a i??* _^ (Ri = 0.5, m = 38 and B = 4) 13 2.4 Basic flow for r* > 1 and B > (Ri = 0.5, m = 38 and B = 2433) 14 2.5 (a) plots the Bingham number versus the rate of viscosity m with Ri fixed at 0.5, (b) plots the Bingham number versus the radius Ri with m fixed at 20 of the Newtonian region 16 3.1 Exact solution for the axy-symmetric basic flow B = 100, m = 38 and Ri = 0.5 36 3.2 Numerical solution for the axy-symmetric basic flow B = 100, m = 38 and Ri = 0.5. 37 3.3 I? norm of the error 38 3.4 Asymmetric basic flow B = 100, m = 38, Ri = 0.5 and (xc, yc) = (0.1,0) 39 3.5 Asymmetric basic flow B = 100, m = 38, Ri = 0.5 and (xc, yc) = (0.2,0) 39 3.6 Asymmetric basic flow B = 100, m = 38, Ri = 0.5 and (xc, yj) = (0.3,0) 40 3.7 (a) Asymmetric flow with B = 4, m — 10, Ri = 0.5 and (xc,yc) — (0.2,0); no unyielded region, (b) Asymmetric flow with B = 100, m = 38, Ri = 0.5 and (xc,yc) = (0.2,0); unyielded region surrounding the Newtonian region, (c) Asymmetric flow with B = 325, m — 5, Ri — 0.5 and (xc,yc) — (0.2,0); yield stress not exceeded 42 3.8 (a) Asymmetric flow with B = 50, m = 38, Ri = 0.5 and (xc,yc) = (0.2,0) , (b) Asymmetric flow with B = 150, m = 38, Ri = 0.5 and (xc,yc) = (0.2,0) 43 3.9 (a) Asymmetric flow with B = 100, m = 10, Ri = 0.5 and (xc,yc) = (0.2,0), (b) Asymmetric flow with B = 100, m = 60, Ri = 0.5 and (xc, yc) = (0.2,0) 45 3.10 (a) Asymmetric flow with B = 100, m = 38, Ri = 0.4 and (xc,yc) = (0.2,0) , (b) Asymmetric flow with B = 100, m = 38, Rt = 0.512 and (xc,yc) = (0.2,0) 47 4.1 Geometry of the perturbed flow 53 4.2 The domain Q,2,c 55 4.3 Coordinate system (r1,9') 57 vii Acknowledgements The author wants to thank to Dr. Ian Frigaard for his invaluable guidance and enormous patience, without him this thesis wouldn't be completed. To Tatiana, for her courage, strength and everything she gave me. To my family and La Morena for their moral and spiritual support. To the Consejo Nacional de Ciencia y Tecnologia (CONACYT) for the financial support. To Dr. Cherif Nouar for all the helpful advices and finally to the readers of this thesis, Dr. Michael Ward and Dr. Brian Seymour. viii Chapter 1 Introduction 1.1 Problem description This thesis concerns the nonlinear stability of a multi-layer flow of two iso-density immiscible fluids in a pipe. Both fluids are driven in the same direction by the action of a pressure gradient, so that we consider a Poiseuille flow. Our basic flow is an axisymmetric flow in which a Bingham fluid surrounds a Newtonian fluid. The flow rates and rheologies are such that the Bingham fluid is unyielded at the interface. This configuration is illustrated in Figure 2.1. A Bingham fluid is a fluid which will not deform under the action of a sufficiently small applied shear stress. Instead, the shear stress must exceed a critical value, known as the yield stress, Ty, in order for the fluid to deform and flow. Bingham fluids are one of the simplest models of a visco-plastic fluid. Some examples of visco-plastic fluids are toothpaste, drilling mud, nuclear fuel slurries, mayonnaise, ketchup and blood. 1.1.1 Applications A principal interest in studying instabilities of multi-layer flows of immiscible fluids is because of their industrial relevance. Two example applications of multi-layer flows are co-extrusion processes and lubricated pipelining. • Co-extrusion processes: Co-extrusion processes are extensively used in the production of bilayer and multi-layer films, sheets and pipes in the plastics industry, the concentric coating of two or more polymers and the production of conjugate fibers in the fiber indus-try. The popularity of multi-layer extrusion has also grown due to its ability to produce 1 Chapter 1. Introduction composite products with tailor-made properties. A co-extrusion operation consists essen-tially of combining several melt streams in a feedblock. The combined melt stream then flows to the die where the layers take their final dimensions. In all co-extrusion opera-tions the rate of production is limited by flow instabilities and especially instabilities at the interfaces between adjacent layers. These instabilities can result in high scrap rates, or in final products with substandard mechanical, optical or barrier properties. • Lubricated pipelining: There is a strong tendency for two immiscible (iso-density) fluids to arrange themselves so that the low-viscosity constituent is in the region of high shear, i.e. typically closest to the wall in a duct flow. Therefore, we can imagine that it may be possible to introduce a beneficial effect in any flow of a very viscous liquid by introducing small amounts of a less viscous lubricating fluid. This is indeed the case in many applications. This process, called lubricated pipelining, is used for example where a heavy crude oil is transported along pipelines with the addition of a small amount of water. In order to control and use any of these multi-layer flows in an industrial process, we must first determine where the two fluids will be. Thus, the dual problems of the stability of interfaces and basic flow configurations are ultimately of great interest. Common to both the above appli-cations is the need to avoid instabilities at the interface between the fluids. These instabilities and their onset often severely limit the rates of production or transport. 1.2 Background literature 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 the latter are of more direct relevance here, the interfacial instability is similar in both flows. 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 [5]. The earliest theoretical investigations of instability due to viscosity stratification is probably 2 Chapter 1. Introduction the classical study of Yih [13]. 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, [6], 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 [12] consider a plane Poiseuille flow, in which 2 fluids flow together down a channel. They extend Yih's results [13] 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. The axisymmetric flow that we consider corresponds, in the Newtonian fluid case, to what is known as perfect core annular flow (PCAF). P C A F has been the subject of fairly extensive study and is reviewed thoroughly in [9]. Here we only point out only a few key details. Hickox [4] studied the linear stability of P C A F using Yih's long wavelength method. The most important conclusion of his work was that no situations were encountered where the primary flow was stable to asymmetric and axisymmetric disturbances, simultaneously. Thus, parallel flows and countercurrent flows are always unstable in the sense that a small perturbation to the flow will initially grow with time at an exponential rate, regardless of the size of the Reynolds number. The most important single cause of instability Hickox found is the difference in viscosity between the fluid regions. Density variation and surface tension have less pronounced effects. However, Hickox's study is confined to the case in which the more viscous fluid abuts the pipe wall. The lubricating case of the less viscous fluid at the pipe wall has been considered by Joseph et ah, [8], who show that stability is possible for a restricted range of parameters. To summarise these results, (and those of many other authors), the commonly accepted method of how to stabilise P C A F , (see [9]), is that it is necessary to place the less viscous fluid at the wall, in order to stabilise against long wavelength disturbances, and to introduce sufficient surface tension in order to suppress short wavelength disturbances. In all cases we note that it is an interfacial mode that is responsible for the instability. 3 Chapter 1. Introduction All the above work has been done for Newtonian fluids. Pinarbasi and Liakopoulos [11] per-formed a linear stability analysis of the interface between two inelastic non-Newtonian fluids in a straight channel driven by a pressure gradient. They considered both Bingham-like fluids and Carreau-Yasuda fluids. They have shown that replacing the bottom Newtonian fluid with a visco-plastic fluid has a stabilizing effect on the fluid-fluid interface for intermediate and large wave numbers. However, Pinarbasi and Liakopoulos do not model a true Bingham fluid with a yield stress. Instead they choose a form of regularised viscosity model in which the effective viscosity attains a Newtonian limit of high viscosity at zero shear-rate. Thus, the results in [11] are best interpreted as effects of having a nonlinearly viscous fluid. The results in [11] are qualitatively similar to those found for the corresponding Newtonian fluid flows. In this thesis we will use a (true) visco-plastic fluid in place of a purely viscous lubricating fluid. Frigaard [2] used a visco-plastic fluid as the lubricating fluid in the study of a three-layer plane Poiseuille flow. The fluids are immiscible, but with no stabilising surface tension. He has shown that linear interfacial instabilities, that would otherwise occur between fluid layers, can be wholly eliminated. This happens when an unyielded plug region is maintained adjacent to the interface. These flows first become linearly unstable to shear instabilities that are analogous to these found for single fluid flows. In this thesis we consider a parallel multi-layer flow of a Newtonian fluid and a Bingham fluid, using the latter as a lubricant. Following the method used in Frigaard [2] and adapting the energy methods in [10], we show that the flow is non-linearly stable for sufficiently small Reynolds number, assuming that the shear stress perturbation can be bounded. This latter assumption is realistic for a well controlled process. 1.3 Outline of the thesis In chapter two we introduce the general model and explain the axisymmetric basic flow. In chapter three we calculate asymmetric basic flow, which are necessary for our stability analysis. We do this using an Augmented Lagrangian method combined with a Finite Element Method. 4 Chapter 1. Introduction We give the outline of the method first and afterwards present several numerical results. In chapter four we consider nonlinear stability via the Reynolds-Orr equation. Although nonlinear our stability bounds are conditional. The conditions we impose correspond to the a priori assumption that the plug region surrounding the Newtonian fluid remains intact for which it is necessary to bound the perturbation of the deviatoric stress. We also allow that the plug and Newtonian fluid region can move. This introduces some complications. We derive three-dimensional nonlinear stability bounds on the decay of the velocity and also a bound for the motion of the plug away from axisymmetry. Chapter 5 concludes with a discussion, conclusions and future work. 5 Chapter 2 Viscoplastically lubricated flows We consider a multi-layer flow of 2 fluids in an infinite pipe, under the action of a pressure gradient; see Figure 2.1 (a). We will focus on the case where the pipe cross section is separated into two distinct fluid domains, with fluid 2 providing a lubricating layer for fluid 1, i.e. the cross-sectional domain occupied by fluid 1 is completely surrounded by fluid 2, which abuts the wall of the pipe. It will be assumed throughout that fluid 1 is a Newtonian fluid and fluid 2 is a Bingham fluid. For simplicity, we consider an arbitrary but finite length L of the pipe and later will consider nonlinear stability of the basic flows to perturbations that are L-periodic with respect to the axial direction. We denote the fluid domain by 0 and the 2 individual fluid domains by Q,\ and 0.2, respectively; see Figure 2.1 (b). Fluid 1 has viscosity / J W and fluid 2 is characterised Theologically by its yield stress r^eld and plastic viscosity /}'2l. Its assumed that both fluids have the same density p. The total flow rate along the pipe is Q and the pipe radius is R, thus defining the mean axial velocity: Uo = Q/TTR 2. The (modified) pressure is denoted p(x,i) and u(x,i) is the velocity. Later, we will use the index k to denote the different fluids. The particular flow that we will consider the instability of is an axisymmetric multi-layered Poiseuille flow. In this configuration there is circular column of Newtonian fluid surrounded concentrically by an annulus of Bingham fluid. For appropriately chosen fluid properties (see §2.2) the Bingham fluid does not yield at the interface and the velocity profile is as shown schematically in Figure 2.1 (c). 6 Chapter 2. Viscoplastically lubricated flows a) Figure 2.1: (a) Newtonian fluid surrounded by a Bingham fluid flowing in a pipe , (b) cross section of the pipe fii corresponds to the Newtonian fluid, Cl2 corresponds to the Bingham fluid, (c) Desired velocity profile (Newtonian fluid surrounded by an unyielded plug region of Bingham fluid) 2.1 Model derivation 2.1.1 Equations of motion The following equations of motion are satisfied within each fluid: dp | dfgW d£i d£j dxj For the Newtonian fluid, the deviatoric stress, is defined by: (2-1) (2.2) fW(«) = A[1]%-(«). (2.3) where (2.4) 7 Chapter 2. Viscoplastically lubricated flows The rate of strain and deviatoric stress second invariants, j(u) and f^(u) k = 1,2 respectively, are defined by: 7(w) 1 1/2 f i 2 i («) 5 E l f M l 2 1/2 (2.5) The constitutive laws for the Bingham fluid use the Von Mises yield condition. They are given by: T^(U) < f y , 4> f l 2 l ( u ) > f y . t-PluA -7 ( « . 7(1*) = 0 7 y ( « ) (2.6) (2.7) At the walls of the pipe, the usual no-slip conditions are satisfied by the velocity, u = 0. In the axial direction there is either no variation in the velocity field, (basic flows), or the velocity field is periodic (perturbed flow). We consider fluids which are either immiscible, but with insignificant surface tension, or that are miscible, but that the timescales are too short for significant molecular diffusion to take place. Under these conditions both stress (traction) and velocity vectors are continuous at the interface between fluids. 2.1.2 Dimensionless Equations of motion We non-dimensionalize the equations of motion (2.1)-(2.7) with the following scaling: P X 4. X = — , t R tUo R ' u u = P= ~ The dimensionless equations of motion are: dui dv,i dt 1 dxj + 1 dnM ^U0 dp dxi ' i t e P l Ox* dxj 0. (2.8) (2.9) (2.10) where p(x, t) denotes the modified pressure, u(x, t) is the velocity and (x, t) is the deviatoric stress. The dimensionless form of the constitutive laws for the Bingham fluid are 7(u) = 0 <^ > T^\U)<B, 1 + T B 7 ( « ) . iijiu) T®(U) > B. (2-11) (2.12) Chapter 2. Viscoplastically lubricated flows and for the Newtonian fluid: [i] 7y = m-iij The dimensionless parameters appearing above are: B = m = [1] (2.13) which are the Reynolds numbers, the Bingham number and the viscosity ratio. At times we will also make use of iJe'1', which is the Reynolds number for the Newtonian fluid: m [i] The rate of strain and deviatoric stress second invariants, j(u) and r(it), are: -| 1/2 1/2 (2.14) Cyl ind r i ca l coordinates In cylindrical coordinates the equations of motion for the Newtonian fluid are: du du v du v2 du dt or r oQ r dz dv dv v dv uv dv dt or r dv r oz dw dw v dw dw m + u ^ + r m + w J z ' 1 d (ru) 1 dv dw r dr r d8 dz 1 fl d f du\ 1 d2u d2u u 2 dv ReW \rdr\ dr J r2 d62 dz2 r2 r2 dB dp dr1 1 fl d ( dv\ 1 d2v d2v v . 1 du (2.15) + + - ^ + ReW \rdr\ dr J r2 d62 dz2 r2 r2 d6 _ldp rdd' (2.16) 1 fl d f dw\ d2w\ dp ReW \rdr \ dr J r2 d92 dz2 J dz' ^ - 0 (2.17) (2.18) 9 Chapter 2. Viscoplastically lubricated flows The dimensionless equations of motion in cylindrical coordinates for the Bingham fluid are du du vdu_v*_ du _ 1 / l 9 / > [ 2 | \ I 9 r { drf} r g \ ~Bt+Udr + rde r+Wdz ReW\rdr\rTrr) + r 88 + dz r dt dr r e r dz ® dp dr' dv dv vdv uv dv 1 ( I d / 2 [ 2]\ , 1 drf} drf} (2.19) dt^ dr^rdd" r ' itei2! ^ r 2 a r V r 86 dz ldp rde' dw . dw . vdw . dw I ( I d , larfj drH]\ dp (2.20) d w _ 1 / }_&_( [2}\ + U~Br~+r~dT+W~dz' ~ \r^dr'\rTzr) dt d ^ r d9 ' ^ dz ReW ^ r 2 a r V 2 r / ' r de ' dz J dz' (2.21) l_d_^+l_dv dw = Q The components of the rate of strain in cylindrical coordinates are: Irr = 2 ^ ) , (2.23) - - <>•»> 2.2 Axisymmetric Basic flows Although will be necessary to consider other basic flows in the course of our analysis, (see §3.1), the key basic flow that we consider is that where the Newtonian fluid is completely surrounded by unyielded Bingham fluid, i.e. a region where the Bingham fluid behaves as a solid. In particular, assume that the Newtonian fluid lies within a circular tube surrounded by an annulus of unyielded Bingham fluid, centered at (x,y) = (xc,yc) = (0,0). This the ideal case. 10 Chapter 2. Viscoplastically lubricated flows The basic flow is calculated assuming that the velocity of the fluids in the radial and angular di-rections is zero and the velocity on the z-direction has r-dependence only, i.e. u = (0,0, W(r)). The Newtonian fluid occupies r G [0, Ri] and the Bingham fluid occupies r G (Ri, 1]. Substi-tuting u = (0,0, W(r)) in the equations of motion we obtain for the Newtonian fluid: 0 = dr' 0 = 0 - -^E. m 1 9 ( y d W l \ dz ReW r dr\ dr J ^ 9 (T9Wl \ = v d p R e [ 2 ] -d r \ d r ) dz m where f\ = —1| and W'(0) = 0, because of symmetry. Without loss of generality we assume / i > 0, so that the pressure gradient drives the fluids in the positive z — direction. Straightforwardly: dW _ r dr h 2 r 2 W(r) = -h—+Cx for 0 < r < Rt. For the Bingham fluid we have 0 0 dp dr'' dp dp l (l d . . dz Re&\ \r dr P r => rrz = -/2-where / 2 = — ff-Re'2' > 0. Thus, the absolute value of the stress \TTZ\ increases linearly with r. It is clear that there are 3 different cases to consider: 1. \rrz\ > B Vr € (r*, 1] and | T „ | < B Vr G (Ri,r*)\ (Case 1). 2. | r r 2 | < B Vr G (-R,, 1]; (Case 2). 3. | T « | > 5 Vr G (-Ri, 1]; (Case 3). 11 Chapter 2. Viscoplastically lubricated Bows We consider the first of these possibilities, since it it this which leads to the type of flow shown schemetically in Figure 2.1 (c). We know that if | r r z | > B B 7~T2. 1 + W2 thus W'-B = -h-We denote the yield surface by r = r*, where rTZ = -B, and hence W'(r*) = 0. Then * 2B B = h r = h (2.29) (2.30) Therefore W = -f2- + f2— + C2 Using the no-slip condition at the wall, W(l) — 0, so „2 r*r j2 r* g W2(r) = - / 2 i r + / 2 ^ + ^ _ / 2 ^ = ^ _ ( ( i _ r ^ _ ( r _ r * ) 2 ) ( 2 . 3 1 Finally, using continuity of the velocity at the interface to evaluate the constant C i , and sub-stituting fi = & = p ^ , the velocity profile is: W(r) = { B ' 1 (R\ - r 2 ) + (1 - r*)2] 0<r<Ri, 2r*lm 2r B_ ^ 2? < r <r*, (2.32) [(1 - r*Y - (r - r*Y\ r* < r < 1. In Figure 2.2 we show an example velocity profile for the ideal case. As we can see, we have unyielded Bingham fluid surrounding the Newtonian fluid. The other two are treated similarly. The velocity profiles are: Case 2: see example velocity profile in Figure 2.3. W(r) = { B r 1 [~(Rf - r 2 ) + (1 - r*f - (Ri - r*)2] 0 < r < ifc, 2r*Lm ^ [ ( l - r * ) 2 - ( r - r * ) 2 ] (2.33) Ri<r<l. 12 Chapter 2. Viscoplastically lubricated Bows 2.5, 2h Figure 2.2: Basic flow for Rt < r* < 1 (Ri = 0.5, m — 38 and B = 100). 13 Chapter 2. Viscoplastically lubricated flows Figure 2.4: Basic flow for r* > 1 and B > ^ (ifc = 0.5, m = 38 and B = 2433). Case 3: see example velocity profile in Figure 2.4. W(r) = { ±[±(R*-r2)] 0<r<Ri, 2r* m (2.34) 0 Ri<r<\. 2.2.1 Determing r* Keeping in mind that our equations have been made dimensionless using the mean velocity, the velocity profile W(r) must satisfy 1 r1 - = J rW(r) dr. (2.35) We now determine r* = r*(Ri, B, m) by satisfying the flow rate constraint. Considering case 1, where we have a velocity profile of type (2.32), after some algebra we derive: 1 t 1 W ^ A B rRi , (l~r*)2(3 + 2r* + (r*) 2) 1 2 = J0 r W ( r ) ^ = 2 ^ [ 4 t + 12 ]' ( 2' 3 6 ) 14 Chapter 2. Viscoplastically lubricated flows which corresponds to find the roots of the following quartic polynomial 0 = (r*)4_4r*(i+ ! ) + 3(i + M). (2.37) Jo m We note that this simplifies to Buckingham's equation for the flow of a single Bingham fluid in a pipe, as Ri —> 0. For Case 2 we have r* < Ri (i.e. no unyielded region). Our basic flow is given by (2.33) and our polynomial is 0 = i2?(4r* - 3Ri) - 4r*(l + | ) + 3(1 + ^ ) , (i.e. continuous with (2.37) as r* —> Ri), which simplifies to the linear equation 0 = r * ( 4 i i ? - 4 - ^ ) + 3 ( l - 3 - i t f ) , D vn Note that since r* < Ri, the boundary of this parameter space is attained when g < VZRj i *? (£ + l ) - 4 i ? i + 3' For Case 3 we have r* > 1, (i.e. the yield stress is not exceeded), the basic flow is (2.34) and our polynomial is simply 0 = 4r*m + BRf (i.e. (2.37) is continuous with it as r* —> 1). Due to the fact that r* > 1, the boundary of this parameter regime is B > | J . (2-39) In Figure 2.5 we explore the dependence of the type of solution on the three parameters Ri, B and m. We note from (2.38) that the boundary between fully yielded and partially yielded regimes, satisfies the following asymptotic limits Am B < ^ as m 0 (2.40) Ri B ± i T T O 8 8 ( 2 ' 4 1 ) B < 4Ri as Rt -> 0 (2.42) 15 Chapter 2. Viscoplastically lubricated Hows 70 m 2001 1 1 r 0 ' 1 "i 1 1 1 1 1 1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Ri Figure 2.5: (a) plots the Bingham number versus the rate of viscosity m with Ri fixed at 0.5, (b) plots the Bingham number versus the radius Ri with m fixed at 20 of the Newtonian region 16 Chapter 2. Viscoplastically lubricated flows Figure 2.5 (a) shows typical variations with B and m, for fixed Ri = 0.5. Figure 2.5 (b) shows typical variations with B and Ri for fixed m = 20. From (2.38)-(2.43) and Figure 2.5 we note that the two constraints (2.38) and (2.39) coincide both as m —* 0 and as Ri —> 1, i.e. limiting the regime where our later results can be valid (since we cannot find the desired basic flow). We also note that the right hand side of (2.38) is maximized V m > 0, as m —> 00, by (2.41), and V Ri € [0,1] by its value at Ri = — ^ — r i i - e - if ( m + ^ ^ 3 B > ; > 4m as m —> 00 ( & + ! ) * - ! we cannot have fully yielded Bingham fluid regions. Intuitively, we expect to require large B and large m, (i.e. through small /t'2'), in order to lubricate. It is clear that this requires some attention to the interface position Ri. 17 Chapter 3 Numerical solution of the asymmetric basic flow The aim of the thesis is to establish nonlinear stability bounds for axisymmetric basic flows of the type illustrated in Figure 2.1c, (see also Figure 2.2). The conventional way of deriving such bounds, (see e.g. [7]), is to consider the evolution of a perturbation about a particular basic flow, such as (2.32), via the Reynolds-Orr energy equation. In general, this is the approach we follow in §4. However, here we have a multi-fluid flow and we must question whether or not the interfacial position remains unchanged when the velocity is perturbed. In [2] a linear perturbation was considered, of the same type of flow as we consider. For a linear perturbation, the linearised momentum equations for the motion of the plug include only linear combinations of the normal modes. Since these equations are integrated over the length of the plug, the net contribution to the plug motion is zero and the plug regions are not perturbed. Unfortunately, for a nonlinear perturbation the same will not be true and we must expect that the normal situation will be that the plug region, (i.e. the interface, which encloses the Newtonian region), moves. However, we wish to preserve the unyielded fluid at the interface, to ensure the stability of the interface. Thus, we later consider a type of conditional instability, in which we bound the perturbation of the stresses \rij(U + u) — Tij(U))\. For such bounded flows we can expect that an annular region about f2i remains unyielded, i.e. even when the Newtonian region moves. For the above type of perturbation, since the Bingham fluid does not yield at the interface, the shape of the interface is unchanged, i.e. fi\ = Qi(t) remains a circular tube of radius Rf, the interface moves but does not deform! 18 Chapter 3. Numerical solution of the asymmetric basic how For these flows, it is natural to perturb about a different basic flow. We denote the centre of Q.\(t) by (x,y) = {xc(t),yc{t)) and seek an axial solution u — (0,0, W(x,y; xc,yc)). This problem is formulated in §3.1. Note we have treated the (axisymmetric) case (xc,yc) = (0,0) in §2.2. Numerical solution is necessary for these asymmetric flows and this is the focus of this chapter. It is clearly important for this particular problem to accurately represent yielded and unyielded regions of the flow. The only method known to do this for flows of a single Bingham fluid is the augmented Lagrangian method, see [1]. Here we generalize this approach to our multi-fluid problem and implement it using the Finite Element Method (FEM). The augmented Lagrangian method is introduced in §3.2. Application of the method to the problem in §3.1 is covered in §3.2.3, and presentation of example results in §3.3. 3.1 Asymmetric basic flows As explained above, we would like to calculate a basic flow u = (0,0, W(x, y; xc, yc)), where the notation means that fii is a cylinder of radius Ri, centered at (x,y) = (xc,yc). In this section we consider cartesian coordinates (x,y,z), and align the z-axis with the axis of the pipe. We assume that \xc\ < 1 — Ri and \yc\ < 1 — Ri so that fii remains embedded in f^ 2-The simplified steady momentum equations are dp _J_ dz + itePl dx dy , (x,y)£Qk k = 1,2 For the Newtonian fluid this simplifies to drzxW drzvW ,d2W d2W + dy dx so that (3.1) is Poisson's equation. For the Bingham fluid, the constitutive laws simplify to: |VW| = 0 <=• [(rz2Jf + ( r ^ < B , 1 + ^ 1 VW [(r] 2J) 2 + ( r | 2 ] ) 2 ]5>5 . | W | J (3.1) (3.2) (3.3) 19 Chapter 3. Numerical solution of the asymmetric basic flow so that we can consider (3.1) as a type of nonlinear Poisson's equation with singular behavior as \VW\ -»0. No-slip conditions are satisfied at the wall: W(x,y) = 0 at x2 + y2 = 1. At the interface (x — xc)2 + (y - yc)2 = Rf, W(x,y,xc,yc) is continuous, as is the traction vector: T£X nx + r]y ny for k = 1,2. The problem described above is one in a wide class of problems studied in [3] using variational methods. In [3] it is shown that for each interface (i.e. here each pair (xc,yc)), there exists a unique weak solution to this problem for given constant ^ . Furthermore, it is also shown that the flow rate through the duct increases monotonically with | | | | . Thus, a unique solution can be found that satisfies the flow rate constraint: 7 r = / W(x,y,xc,yc) dxdy. (3.4) Jx2+y2<l In practice this solution must be found numerically, by iterating with respect to | ^ |. 3.1.1 Variational form of the problem In order to solve numerically the basic asymmetric flows two methods can be used. (i) We can regularize the effective viscosity. In this method, the singularity in the effective viscosity is replaced by a Newtonian limit at zero shear rate, e.g. one option is: B 1 + VW (3.5) | V W | + eJ with e -C 1. It can be seen that using this method we are not modelling a true Bingham fluid, at low shear rates. In its favour this method can be implemented in most stan-dard codes that allow a nonlinear viscosity. Additionally, since the effective viscosity is differentiable, faster iterative methods may be used as a part of the numerical solution. One major drawback is that one needs to examine convergence with respect to e and then decide what value of e is "small enough". Secondly, it must be remarked that there are 20 Chapter 3. Numerical solution of the asymmetric basic flow-many choices for a regularization such as (3.5), and it is unclear which should be used. Finally, in a problem such as ours, we specifically wish to determine the yield surface, or at least to know definitively that the fluid is unyielded. Therefore, we reject this method. (ii) The augmented Lagrangian method is a method that allows the zero shear rate to be eval-uated exactly, hence computing the unyielded region correctly. This is an overwhelming advantage for our problem and therefore we adopt the augmented Lagrangian method. We give an outline of the method in §3.2. To apply this method to our problem we need to express it in a variational form. Recall that the simplified momentum equations for the asymetric basic flow, u = (0,0, W(x, y; xc, yc)), are Bz iteP] Br [fcl Br [fcl U'zx _j_ ulzy 1,2, (3.6) Bx ' By Multiply (3.6) by v € HQ(Q) and integrate over the whole domain f fv dx + f (V • T)V dx = 0 (3.7) JQ JQ where r = (Tzx,Tzy) k = 1,2. Assuming without loss of generality that / > 0, so that the pressure gradient drives the fluids in the positive ^-direction, / = — gf-Re'2' = J?e'2'|g2|. Integrating (3.7) by parts we have: fvdx+ ( T • n)v ds — / T • Wv dx, = 0 for k — 1,2. JQ. Jduk JQ, The boundary integrals vanish because v 6 HQ (Q), and of the continuity of the traction over the interface, i.e., [r^ • rij] = 0 for k, j = 1,2. Thus, we are left with: f fvdx- f T • Vu dx, = 0 (3.8) JQ JQ which, on following [3], becomes the following variational inequality: - ( f(u- W) dx + uk f VW- (Vu - VW) dx jQk JQk + Bk f |Vu | - | V W | dec > 0. forfc = l,2. (3.9) JQk 21 Chapter 3. Numerical solution of the asymmetric basic flow-As in [1], this is equivalent to the following minimization problem min J(W), (3.10) where J(W) = [ E ^a(W, W) + Bkj(W)] - L(W), fc=l,2 a(W,W)= f VWVWdx, j(W)= f \VW\dx, for ft = 1,2 L(W) = f fW dx, Hi -m, [12 = 1, 0, B2 = B, Unfortunately, we can not directly apply Newton's method or conjugate gradient to our min-imisation problem (3.10) because of the non-differentiability of jiW). We note that we could also use a regularisation method here, by letting with e <S 1, but this is more or less equivalent to solving our original problem by substituting a regularised viscosity. Therefore, we will use the augmented Lagrangian method, (see [1]), which is a method for solving the non-differentiable minimization (3.10), rather than a perturbed form of the classical problem (3.6). 3.2 Augmented Lagrangian Method In this section we define what an augmented Lagrangian is, we state the type of problems that lead us to this representation and we give an algorithm to solve the minimization problem (3.10). 22 Chapter 3. Numerical solution of the asymmetric basic flow 3.2.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 d £ R . We associate the following quadratic functional to these parameters J(v) = ±(Av,v)-(b,v) (3.11) where in (3.11), (•, •) is the canonical Euclidian inner product in R ^ . Let B be a linear mapping from R ^ into HM, this can be identify by a M x N matrix. Now we consider the minimisation problem J(u) < J(v) G Ker B (3.12) u G KerB. (3.13) It is known that (3.12) has a unique solution, but the constraint (3.13) may make the solution difficult to compute. Introducing a Lagrange multiplier p to (3.12), in order to get an unconstrained problem, we have min {J(v) + (p,Bv)} (3.14) v G RN 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,Bv) (3.15) recall that (u,p) will be a saddle-point of L if L(u,q) < L(u,p) < L(v,p) W G R ^ , q e R M (3.16) 23 Chapter 3. Numerical solution of the asymmetric basic flow thus we have (3.17) It can be shown that L admits at least one saddle-point (u,p) on KN x R , where u is the solution of (3.12) and is common to all the saddle-points of L on R ^ x R M . We define the augmented Lagrangian Lr , for r > 0, by and note that a saddle-point of L(v,q) for which u € KerB, will also be a saddle-point of Lr(v,q). However, (3.18) may be better conditioned to solve numerically. 3.2.2 Algori thms using Augmented Lagrangian Methods We now present an algorithm for finding the saddle-points of an augmented Lagrangian func-tional Lr. We generalize the above problem to include our asymmetric flow problem (3.10). Consider as our generalized minimization problem, where • V and H are normed vector spaces (real for simplicity). • B e L(V, H). • F,G are functions which are convex, proper, and lower semi-continuous (l.s.c) on H and V, respectively. The following algorithm is due to [1], who we follow closely in the sequel. Lr(v,q) = J(v) + (q,Bv) + -\Bv 2 = L(v,q)+V-\Bv (3.18) mm{F(Bv) + G(v)} (3.19) veV 24 Chapter 3. Numerical solution of the asymmetric basic How The essential idea of the algorithm is based on the fact that there is an equivalence between (3.19) and Min {F(q) + G(v)} (3.20) (v,q)€W with W = {{v, q) G V x H, Bv - q = 0}. Introducing the variable q, linked to v through the linear equality equation Bv = q we will utilize a Lagrange multiplier and reduce (3.20) to a saddle-point problem. We will assume that the spaces V and H are Hilbert spaces; we denote by (•, •) the inner product on H, and by | • | the associated norm. We define, for v G V, q G H, n G H the Lagrangian L(v, q, n) = F(q) + G(v) + (/z, Bv - q) and then for r > 0, the augmented Lagrangian is defined as L(v, q, n)r = L{v, q, /x) + ^\Bv - q\2. (3.21) 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 (3.21). A L G 2 (p°, A1) G H x H arbitrarily specified; (3.22) with ( p n - 1 , A n ) known, determine succesively u n , p n , \ n + l by G(v)-G(un) + (Xn,B(v-un)) +r(Bun - p n _ 1 , B ( u - «")) > 0 Vv G V, (3.23) un G V 25 Chapter 3. Numerical solution of the asymmetric basic flow F(q) - F(pn) - (A n, q - pn) + r(pn - Bun, q - pn)) > 0 \/q G H, (3.24) pn GH Xn+1 = ) n + p n ( B u n _ pny (3 2 5 ) Convergence of A L G 2 We will follow the convergence proof shown in [1]. In order to get convergence of A L G 2 (3.22)-(3.25) assumptions over F and G are necessary, (i) B is injective and ImB is closed in H . (3.26) (ii) lim . . = oo M \q\ - oo (3.27) (iii) Taking into account the properties of G, (3.26) and (3.27) imply that (3.19) admits a unique solution u. In the following text we will write p = Bu. (iv) Assume that F = FQ + F\, where F\ is convex, proper and l.s.c. on H, and where Fo is convex, 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 5M • [0,2M] —> strictly increasing, with SQ = 0, such that for all p,q G H, \p\ < M, \q\ < M, we have: (Fi(q)-Fa0(p),q-p)>6M(\q-p\), where, in (3.28), FQ denotes the gradient of Fo. 26 (3.28) Chapter 3. Numerical solution of the asymmetric basic How Theorem We assume that Lr admits a saddle-point (u,p, A) on V x H x H. Under the assumption on B, F, G (3.26)-(3.28) and if pn satisfies 0 < Pn = P < (^Y^)T, (3.29) we have for A L G 2 the following convergence results un —• u strongly in V, (3.30) pn —> p strongly in H, (3.31) A n + l _ Xn _^ o strongiy i n H^ (3.32) A n is bounded in H. (3.33) Moreover, if Xn is a (weak) cluster point of Xn in H, (u, p, A) is a saddle-point of Lr on V x H x H. Proof (following [1]): Let u~n — u n - u , p~n = pn-p, A _ n = A n - A, and (u, p, A) a saddle-point of Lr on V x H x H, we have G(v)-G(u) + {X,B(v-u)) +r(Bu - p, B(v -u))>0 Vv G V, (3.34) (F"0{p),q-p) + F1(q)-F1(p) -(X,q-p) + r{p -Bu,q-p)) > 0 G H, 27 (3.35) Chapter 3. Numerical solution of the asymmetric basic How X = X + pn(Bu-p) Vn. (3.36) from (3.22)-(3.25) we have G(u)-G(un) + (Xn,B(v-un)) +r(Bun - pn~\B(v - un)) > 0 Vu € V, (3.37) (F^pn),q-pn) + F1(q)-F1(pn) -(Xn,q-pn) + r{pn - Bun,q-pn)) > 0 V? € H, (3.38) An+1 = A n + p n ( B u n _ pnj y n (3 39) Setting v = un in (3.34), g = in (3.35) then v = u in (3.37) and <? = p in (3.38); hence by addition r\Bu-n\2-r(p-n-\Bu-n) + (X-n,Bu-n) <0 (3.40) (F6V) - *"(p),Pn - P) + HP"71!2 - r(p-n, Bu'n) -(X~n,p-n) <0 (3.41) Adding (3.40),(3.41) and regrouping the terms, we obtain (F6V) - F^(p),pn -p)+ r\Bu~n - p~n\2 + (A"", Bu~n - p~n) +r(p-n - , Bu~n) < 0 (3.42) Also subtracting (3.36) from (3.39) and taking scalar square in H, we have | A - " | 2 - | A - r a + 1 | 2 = -2p(Bu-n - p~n, X~n) - p2\Bu-n - p-n\2 (3.43) It then follows from (3.42) and (3.43) that | A - » | 2 _ |A-n+l|2 > 2p{F^pn) - F^(p),pn - p) + p(2r - p)\Bu~n - p-n\2 + 2pr(p-n -p-n-\Bu-n) (3.44) 28 Chapter 3. Numerical solution of the asymmetric basic flow We will try to modify the last term on the right hand side of (3.44); starting from Bu~n = (Bu-n - Bu'71-1) + (Bu'71-1 - p-n~x) + p-71-1 Thus, (Bu-n,p-n -p~n'x) = (Bu~n - Bu-n-\p-n -p~n~l) +(Bu-n~1 - p-n-\p-n - p-71'1) +{P~n~\v~n ~ P~n~l) (3-45) Prom (3.45) and from (p-n-\p-n-p-n~1) = \(\p-n\2 - IP""-1!2 - \p~n-p-71'1?) we deduce that 2pr(Bu-n,p-n - p ~ n _ 1 ) = 2pr(BvTn - BuTn~l,p~n - p~n~l) +2pr(Bu-n-1 - p-n-\p-n - p'n'x) +pr(\p-n\2 - b""- 1 ] 2 - \p'n-p-n~l\2). (3.46) Considering (3.38) on iteration n — 1 instead of n, we have ( W _ 1 ) , < Z - P n ' 1 ) + Fi(q) - F^-1) -(Xn,q-pn-1) + r(pn-1 - Bu^^q-p71-1)) > 0. (3.47) Taking q = pn~x in (3.38) and q = pn in (3.47) we obtain by addition (*o(pn) - ^ ( p n - 1 ) , p n -p11-1) - (X~n - X~n~l,p~n - p~n~l) +r\p-n - p~n-l\2 - r(Bu-n - Bu-n-l,p-n - p"""1) < 0 (3.48) It follows from the monotonicity of FQ, and from (3.48), that r\p~n - P " n _ 1 | 2 - ( A " n - X-71-1 ,p~n - p-71'1) -r(Bu-n - Bu-n-\p-n - p'71'1) < 0 (3.49) 29 Chapter 3. Numerical solution of the asymmetric basic Row We have from (3.25) It then follows from (3.49) and (3.50) that r(Bu-n - Bu-n-\p~n-p-"-1) > r\p~n - p-n~l\2 - p{Bu~n-1 -p-n-\p-n -p-"-1). (3.51) It then follows from (3.46), (3.51) that 2pr(Bu-n,p-n - p-n~l) > pr(\p-n\2 - \p-n~1\2 + \p~n - p-71'^2) + 2p{r - p){Bu-n-1 -p-n-\p-n -p-n~l). (3.52) Finally, combining (3.44), (3.52), we obtain ( | A - f + pr\p-n-l\) - (\\-n+1\2 + pr\p-n\) > pr\p-n - p-"-1]2 + 2p(F0'(pn)-F0-(p),pn-p) + p(2r - p)\Bu~n - p~n\2 + 2p{r - p^Bu-71'1 -p~n-x,p-n -p-™-1). (3.53) Using Cauchy-Schwarz inequality in (3.53) we have ( |A-"| 2 + PT-IP-"-1!) - (\\-n+1\2 + pr\p-n\) > pr\p-n - p-™-1]2 + 2p(F^pn)-F^p),pn-p) + p(2r - p)\Bu~n - p-n\2 - p\r-p\(\Bu-n-1-p-n-1\2 + \p-n-p-n-1\2). (3.54) Assuming p = r we obtain | A - f - | A - " + 1 | 2 > 2 r ( ^ ( p n ) - ^ ( p ) , p " - p ) + r 2 | J B « - " - p - " | 2 . (3.55) Thus, we have lim \Bun - pn\ = 0 (since Bu=p), n -* oo (3.56) 30 Chapter 3. Numerical solution of the asymmetric basic flow lim (F^pn)-F^p),pn-p)=0, n —> oo (3.57) and we have A n bounded in H. We will show now that pn is bounded in Ho- In fact, since the functional F is proper, there exists p£H, such that —oo < F(p < oo). Then putting q = p in (3.24), we obtain F(p)-(Xn,p) + r(pn-Bun,p) > F(pn)-(Xn,pn) +r(pn - BUn,pn). (3.58) Since A" is bounded in H and knowing that pn — Bun tends to zero, we thereby deduce that there exists positive constants Po, Pi, independent of n, such that Po > F{pn) - px\pn\. (3.59) Taking into account assumption (3.27), it follows from (3.59) that pn is bounded, i.e. that there exists M > 0, such that \pn\ < M, Vn. (3.60) Furthermore, for M sufficiently large, we also have |p| < M. (3.61) It then follows from (3.60), (3.61) and from the uniform convexity of Fo over the bounded subsets of H, that we obtain (F0-(pn)-F0-(p),pn-p)>6M(\pn-p\) Vn. (3.62) Thus lim 6M(\pn-p\) = 0-n —> oo (3.63) 31 Chapter 3. Numerical solution of the asymmetric basic Bow It then follows from the properties of 5M that lim \pn - p\ = 0. n -> oo (3.64) and from (3.56) that lim Bun=p (= Bu), strongly in H. n —* oo (3.65) Since the operator B is injective with ImB closed in H, this implies that lim un = u, strongly in V. n —> oo (3.66) Now, if 0 < p < r, taking into account (3.54), we have ( |A- n | 2 + prlp-"" 1! + p(r - p)\Bu-n~l - p~n~l\2) _ ( | A -n+l |2 + ^jp-nj + p(T _ p ) \ B u - n - p~n\2) > > 2p(F()(pn) - F^p),pn -p) + pr\Bu~n - p~n\2 + p2\p~n - p-n~l\2, (3.67) which implies, just like when p = r, (3.30)-(3.33). If p > r we have \r — p\ = p — r; it then follows from (3.54) that the convergence results (3.30)-(3.33) will apply if we have p < PM, where PM{2T - PM) = ~PM{PM - r), a pMr = apM(pM - r). (3.68) By elimination of a we deduce from (3.68) that PM ~ RPM - r 2 = 0 i.e. (since PM > 0) l + v/5 PM = r. 32 Chapter 3. Numerical solution of the asymmetric basic flow 3.2.3 Application of ALG2 to our problem. Introducing the extra constraint VW — q = 0 to (3.10) it is easy to see that (3.10) is equivalent to (3.20) with B = V, G(u) = - f fudx, for A; = 1,2 F (9) = X ^ f c / QQdx + Bk \Q\dx, q = Vu. where the first part of F(q) is Fo and the second is F\, H = L 2 (Q) and V = HQ(Q). Therefore the unconstrained problem with its correspondent Lagrange multiplier is mm{G{u) + F(q) + (A, q - VW)}. ueH^qeL2^) And introducing the Augmented Lagrangian defined for r > 0, we have the minimization problem min{G(u) + F(q) + Uq~ VW, q - VW) + (\,q- VW)}. (3.69) W eH^q£L2(^2) We note that B\ = 0 in Qi. Therefore, in Qi we have a linear problem to solve and would like to avoid iterating with respect A and q. Minimization of (3.69) with respect to u is equivalent to solving 0 = - / - rV2u + rV q + V • X (3.70) in each domain Q^, provided that along interface <9fii the quantity / u(rVu — rq — X) • n ds (3-71) is continuous from both sides. 33 Chapter 3. Numerical solution of the asymmetric basic flow To avoid iterating in Cli, we enforce A = q = 0, and to ensure that (3.70) represents (3.1) in Cl\ we set r = m. Note that the algorithm A L G 2 converges for arbitrary r > 0, provided that p satisfies (3.29). By virtue of our discretisation using F E M , with nodes placed on the interface, u will be contin-uous in (3.71). Secondly note that as n —> oo, un -* u and qn —* q with q = V u . Thus (3.71) becomes / umVu • n ds p m j d 1 = / ~ ~ u ^ ' n + w(rVu — rq) • n ds (3.72) Fluid 2 but from (3.70) it is easy to see that A converges to the stress, so that (3.72) just represents continuity of the stresses at the interface. With the above modifications, A L G 2 for our problem is: • Fix qn, A" minimise with respect to u in the following functional, / ~fu + T T | V U | 2 - (mq + A) • V u dx + -fu + ^\Vu\2 dx to find u n + 1 . • Fix un+1, A n minimise with respect to q over the following functional: I ]-^\q\2+B2\q\-(™^un+1-\n)-q dx J9.o 2 m2 to find qn+1. • Fix qn+1, un+1 maximise with A to give A n + 1 . With the continuity condition (3.72) and continuity of velocity u at the interface dili, we are in fact solving 0 = / + m V 2 u x G Cli 0 = / + m ( V 2 - V - q ) - V - A xeCl2 using the F . E . M . method. 34 Chapter 3. Numerical solution of the asymmetric basic flow For the minimisation with respect to q we solve min{(^)\q\2 + B\q\-M-q}, Q where M = A + r V u . Thus, If \M\ < B q = 0 If | M | > 5 => g = <9M, where # = ^ j ^ ^ and B is the Bingham number. Finally, our update for A is as follows: A " + 1 = \ n + p(qn+1 - V « n + 1 ) , where p = r. 3.3 Example Results Here we present results for our computations of W(x, y; xc, yc) using the augmented Lagrangian method. The algorithm has been implemented in MatLab and we have used FEMLab for the Delaunay triangulation of the domain. For the majority of our results we have used 2000 elements and 1000 nodes. We start in §3.3.1 by showing that the method converges, by running test cases against the axisymmetric analytic solution. For the remainder of the section we explore how the basic flow changes with (xc,yc) and with the process parameters B, m, Ri. In considering changes with (xc, Uc) our primary concern is to demonstrate that the plug region is preserved for sufficiently small departures from symmetry, see §3.3.2. For the later examples, as well as showing the velocity profile, we compute the pressure gradient / = i?et2]||2| ; the flow rate in the Newtonian region, Q\, and the flow rate in the Bingham region, Q2- These flow rates Q\ and Q2 are formally defined by Ql = W(x,y) dxdy ^(x-x c)2+(y-yc) 2<H? Q2 = W(x, y) dxdy jR?<(x-xc)2+(y-yc)2<l 35 Chapter 3. Numerical solution of the asymmetric basic Row Exact solution of the basic flow 1-5> -1 -1 Figure 3.1: Exact solution for the axy-symmetric basic flow B = 100, m — 38 and Ri = 0.5. and we recall that: 7T = / W(x, y) dxdy = Q i + Q 2 -Jx2+y2<l Additionally, we tabulate the minimum distance from dfl\ to the yield surface, where it exists. These parameters are those which would be of interest in attempting to control the process. Although we would require a more comprehensive database of parametric variation, the results we present at least allow one to compute the sign of the derivatives of these quantities, with respect to changes in each of the parameters. 3.3.1 Convergence We first demonstrate convergence of the method. Figure 3.1 and Figure 3.2 show the results of computing the basic flow for parameters Ri = 0.5, ro = 38, B = 100. In Figure 3.1 is shown the exact solution, represented on the same triangular elements as Figure 3.2, which shows the results of the augmented Lagrangian method. Figure 3.3 plots the L2 norm of the error of this example as the mesh size h is decreased. As 36 Chapter 3. Numerical solution of the asymmetric basic Bow Numerical solution of the basic flow 1.5 -1 -1 Figure 3.2: Numerical solution for the axy-symmetric basic flow B = 100, m = 38 and Ri = 0.5. we can see, we have a linear decay of the error norm. This is consistent with the convergence results for the Finite Element Method. The error is 0(h) where h is the maximum over the edges lengths of all triangles. We observe in Figure 3.2 that the unyielded plug region surrounding fii is very evident. Al-though this type of flow might be computed using a viscosity regularization (as discussed earlier), the gradient of W within the unyielded region would be non-zero, (in fact, typically of the order of the regularization parameter used). Here, in our implementation of A L G 2 , see §3.2.3, it is very clear that the gradients are set identically to zero. 3.3.2 Preservation of the plug region. Our later analysis relies on the fact that an axisymmetric basic flow of the type shown schemat-ically in Figure 2.1c, when perturbed to an asymmetric basic flow, still preserves a plug region in an annular region surrounding Or. To demonstrate this numerically, we have fixed the parameters B = 100, m = 38, Ri = 0.5, (see Figure 3.1), and computed the basic flow for 37 Chapter 10" o £ 10" I LLI 10' .-3 10 3 Number of Figure 3.3: L2 norm of the error. (xc,yc) = (0.1,0), (0.2,0) and (0.3,0). The results are shown in Figures 3.4-3.6. Whereas for the smaller perturbations the plug region is preserved, as (xc,yc) becomes large (Figure 3.6) the plug is destroyed on the side closest to the pipe wall. To understand this perturbation process slightly better, consider the axisymmetric flow of §2.2. Having found r* via (2.37), the shear stress is given by rj 2 ' = — ^ T y 2 ' and the pressure gradient is | | = — 2 ^ « P 1 TY2^. The flow is axisymmetric, hence TQZ = 0 everywhere and r = |rj 2 ' | = p-ry2\ Thus, at the interface i.e. the larger the initial annulus of unyielded fluid, the further below the yield stress is the interfacial stress, and the more likely that the plug region is preserved. Now consider small |(x c,y c)| ^ 0. If the plug region is preserved then W(x,y,xc,yc) satisfies Re® dp (3.73) m dz in fii where rc = (x2. + y2)K i-e. it is evident that, (apart from a rotation about (0,0)), the 38 Chapter 3. Numerical solution of the asymmetric basic flow Numerical solution of the basic flow, first perturbation tJS>, -1 -1 Figure 3.4: Asymmetric basic flow B = 100, m = 38, Ri = 0.5 and (xc,yc) = (0.1,0). Numerical solution of the basic flow, second perturbation -1 -1 Figure 3.5: Asymmetric basic flow B = 100, m = 38, Ri = 0.5 and (xc,yc) = (0.2,0). 39 Chapter 3. Numerical solution of the asymmetric basic flow Numerical solution of the basic flow, third perturbation Figure 3.6: Asymmetric basic flow B = 100, m = 38, Ri = 0.5 and (xc,yc) = (0.3,0). basic flow depends on rc. When r c = 0 we have §f (0) = — 2 R ^ and it is reasonable to expect that | | varies continuously with r c . The boundary conditions for W in (3.73) are that W = Wp at the interface, where Wp is the plug velocity. Thus, W satisfies Poisson's equation with constant boundary conditions on fij, for which there is a unique solution. Furthermore, this unique solution is axisymmetric, depending only on r' = [(x — xc)2 + (y — y c ) 2 ] 5 -Therefore, if the plug region is intact, inside the Newtonian domain we have T r ' z = ^2c^z~m~ a n d T9'* = ° where 8' is an azimuthal coordinate centred at (xc,yc). Continuity of stress ensures that TT>Z is continuous at the interface r' — Ri, (but not necessarily Tg/Z). Thus, if g is continuous with rc and if \TQ>Z\ varies continuously from its values at rc = 0, (i.e. from zero), then r' 2' < B at the interface for some finite range of sufficiently small rc. On the other hand, since the no-slip condition must be satisfied, it is quite clear that the plug must 40 Chapter 3. Numerical solution of the asymmetric basic flow Fig. B m Ri (Xc, Vc) / Qi Q2 size of the plug 3.7a 4 10 0.5 (0.2,0) 18.6073 1.1325 2.0091 no plug region 3.7b 100 38 0.5 (0.2,0) 231.7366 0.9764 2.1652 0.0862 3.7c 325 5 0.5 (0.2,0) 644.8610 3.1416 -8.9445E-06 0.3 Table 3.1: Comparison of various asymmetric basic flows. yield at the interface as rc —> 1 — Ri and Oi touches the wall. 3.3.3 Different solution types for asymmetric basic flows Here we present a series of examples to demonstrate that the 3 solution types found for the axisymmetric flow (see §2.2) are also evident when the flow is asymmetric, i.e. just perturbed asymmetrically. In Figure 3.7a we show a velocity profile where the unyielded region does not exist. In Figure 3.7b we present the ideal case where an unyielded Bingham fluid region surrounds the Newtonian fluid region. In Figure 3.7c we show a velocity profile where the yield stress is not exceeded and therefore, due to the no-slip condition, the Bingham fluid is not moving at all. The computed flow parameters are presented in Table 3.1. 3.3.4 Variations with B In Figure 3.8 we fix the parameters Ri = 0.5, m = 38, (xc, yc) = (0.2,0), and vary the Bingham number. The results for B = 50 and B — 150 are shown in Figures 3.8a and 3.8b respectively. These can be compared with our base case in Figure 3.7b, which has the same parameters, but for B = 100. Table (3.2) tabulates the computed parameters. We can see that the pressure gradient decreases or increases significantly when we vary the Bingham number. The flow rates Qi and Qi do not change significantly. We can see that when the Bingham number decreases or increases the unyielded plug region becomes thinner or wider as expected. 41 Chapter 3. Numerical solution of the asymmetric basic flow -t -1 Figure 3.7: (a) Asymmetric flow with B = 4, m = 10,Ri = 0.5 and (xc,yc) = (0.2,0); no unyielded region, (b) Asymmetric flow with B = 100, m = 38, Ri = 0.5 and {xc,yc) = (0.2,0); unyielded region surrounding the Newtonian region, (c) Asymmetric flow with B — 325, m — 5, Ri = 0.5 and (xc,yc) = (0.2,0); yield stress not exceeded. 12 Chapter 3. Numerical solution of the asymmetric basic flow Chapter 3. Numerical solution of the asymmetric basic flow Fig. B m Ri (XcVc) / Qi Qi size of the plug 3.8a 50 38 0.5 (0.2,0) 123.552 0.9547 2.1809 0.0506 3.7b 100 38 0.5 (0.2,0) 231.7366 0.9764 2.1652 0.0862 3.8b 150 38 0.5 (0.2,0) 338.1976 1.0133 2.1283 0.1052 Table 3.2: Variations with Bingham number. Fig. B m Ri (xc, Vc) / Qi Qi size of the plug 3.9a 100 10 0.5 (0.2,0) 229.2594 1.2688 1.8728 0.1052 3.7b 100 38 0.5 (0.2,0) 231.7366 0.9764 2.1652 0.0862 3.9b 100 60 0.5 (0.2,0) 231.9616 0.9373 2.2043 0.1420 Table 3.3: Variations with viscosity ratio, m. 3.3.5 Variations with m In Figure 3.9 we illustrate how the asymmetric basic flow changes with viscosity ratio, m. Parameters are fixed at Ri — 0.5, B = 100, (xc,yc) = (0.2,0) and we present results for m = 10 and m = 60; (again compare with Figure 3.7b, where m = 38). In Table (3.3) we can see that the pressure gradient increases slightly as m increases, as is physically intuitive. Similarly, the Newtonian flow rate drops, since it becomes relatively easier to flow in the Bingham region. The effect on the size of plug region is ambiguous. Physically, as we increase m, if Q\ were fixed, then we would increase the shear stress in the Newtonian region and hence decrease the plug. However as we have seen, increasing m reduces Q\ and therefore the net effect on the shear stress is unpredictable. 3.3.6 Variations with Ri Finally, Figure 3.10 looks at variations with Ri. Parameters are fixed at m = 38, B = 100. The results for Ri — 0.4 and Ri = 0.512 are shown in Figure 3.10a and 3.10b, respectively (c.f. Figure 3.7b again). In Table (3.4) we can see that the pressure gradient decreases with Ri in this particular case. However, this is not conclusive. For the case presented the Newtonian 44 Chapter 3. Numerical solution of the asymmetric basic flow Chapter 3. Numerical solution of the asymmetric basic flow Fig. B m Ri / Qi Q2 size of the plug 3.10a 100 38 0.4 (0.2,0) 232.3296 0.6066 2.535 0.1955 3.7b 100 38 0.5 (0.2,0) 231.7366 0.9764 2.1652 0.0862 3.10b 100 38 0.512 (0.2,0) 226.6466 0.9929 2.1487 0.0806 Table 3.4: Variations with R4. fluid is probably less viscous than the Bingham fluid, (i.e. considering the effective viscosity of the Bingham fluid). Were the situation reversed, the we may have seen an increase in pressure gradient. An increase in Q\ is to be expected as is the decrease in size of plug region, since the entire Bingham region decreases. 46 Chapter 3. Numerical solution of the asymmetric basic flow Chapter 4 The Reynolds-Orr equation and Stability bounds In this chapter we derive the nonlinear stability bounds, that are the aim of the thesis. Although nonlinear our bounds are conditional rather than global. The conditions we impose are the following: (i) The basic flow centered at (xc,yc) must have an intact plug region surrounding f2i(t). (ii) For any basic flow satisfying (i) we bound the perturbation of the deviatoric stresses. In §3.3 we have seen that when the axisymmetric basic flow (2.32) has a plug region surrounding the Newtonian fluid, then this plug region still remains for asymmetric basic flows, provided that |(x c,y c)| is sufficiently small. Thus, condition (i) is reasonable. For (ii), we expect that by taking the bound sufficiently small, the plug region will persists in an annular region surrounding the Newtonian fluid. Thus, for small finite perturbations, the Bingham fluid cannot deform close to the interface with the Newtonian fluid. Since the pipe is considered infinite and the unyielded Bingham fluid cannot touch the walls, it is only possible for the plug to have a linear motion, plus a rotation about an axis parallel to the z-axis. The Newtonian region f2i(i) remains a circular tube, of radius B4, that at time t has axis centred at (xc(t),yc(t)). It is the stability of this flow that we consider below. 48 Chapter 4. The Reynolds-Orr equation and Stability bounds 4.1 Derivation of the Reynolds-Orr equation Regardless of whether or not the plug region remains intact, there exists a unique solution to the steady axial flow problem considered in §3.1. In §3.2.3, we have demonstrated numerically that the plug region persists in the basic flow, for small finite \(xc(t),yc(t))\, i.e. when the solution to the axisymmetric problem permits a plug region surrounding fijv- We therefore consider perturbed velocity and pressure fields of the following form: U + u = ( 0 ,0 , W(x, y; xc(t),yc(t)) + (u,v,w), P + p = P(z; xc(t),yc(t)) + p, where we explicitly state the dependance of the basic flow on (xc,yc), as in §3.1. Thus, our basic flow is a full solution to the steady Navier-Stokes equations, but is not itself steady. Time dependence in (U,P), enters only from the rigid motion of the plug, i.e. changes in (xc(t),yc(t)). In our derivation of stability bounds, we first derive conditions under which the perturbation, | |u| | , decays. We then show that the decay in allows only a finite movement of the plug region, i.e. \(xc(t),yc(t))\ is bounded for all time. By choosing the rheological parameters appropriately and bounding the initial perturbations this verifies our a priori assumption regarding the preservation of the plug around Qi (t). For the basic and perturbed flows, at time t the fluid domains are fii(i) and ^( i ) , a n d the equations of motion in fluid k are: Oi + ^ + uj) — (Ui + Ui) = - ^ P + p ) + ^ ^ i U + u). (4.3) o - ^ T ( ^ + 0 (4-4) (4.5) Note as above, that strictly speaking (4.1)-(4.2) are not the equations of steady motion, since (U, P) are not steady. Secondly, note that by defining (U, P) with relation to (xc(t),yc(t)) the 49 Chapter 4. The Reynolds-Orr equation and Stability bounds fluid domains of the basic and perturbed flow coincide at time t. This is the primary motivation for perturbing around W(x,y;xc,yc), rather than W(x,y;0,0). The point (xc(t),yc(t)) is moving with speed (uc(t), vc(t)), given by: uc(t) = ^ (t), vc(t) = ^ (t). (4.6) This allows us to rewrite (4.3) as follows: di + ^ + U ^ d x -. rdW dW , .. -6i3[dx-cUc+dy-cVc]' (4-7) where 5ij is the Kronecker delta function. Note that — 0. We subtract (4.1) from (4.7), multiply by ttj, and integrate over the individual fluid domains, at time t: f D f . dW dW dW . , j^-Uidx = -^[uj— + —uc + — v c ] W d X +L u ^ - p + R W ] [ T ^ i u + u ) ~ T ^ m ] n ^ ] d s ' (4-8) where we have used the divergence theorem to derive the last term; denoting the outward normal to dQfc. In (4.8), ^ denotes the material derivative. We now sum the 2 fluid domains in (4.8). Note that all the boundary integrals vanish, through a combination of: (i) At the wall u\ = 0, (ii) At the ends of the flow domain, say z = ±L/2, the perturbation is periodic and the outward normals change orientation, (iii) At the interface both u and the traction vectors of both basic and perturbed flows are continuous. 50 Chapter 4. The Reynolds-Orr equation and Stability bounds Thus, finally we have: D_ Dt - f dx dW | dW '3 dxj dxc dW . uc + —— vc]w ^ E ^ ^ T ^ U + « ) - to- (4-9) which is the Reynolds-Orr equation for the kinetic energy of the velocity perturbation u. We analyse this equation further in §4.2, to derive our conditional bounds. Note that in considering [T\^ (U+U) — T^ (U)] at a point (x, t), since the fluid domains coincide for the perturbed and unperturbed flow it suffices to expand the constitutive laws for fluid k. If instead, W(x,y,0,0) had been taken as the basic flow then the fluid present at (x,t) for perturbed and basic flows can be completely different. It then becomes difficult to bound This problem is quite classical for any multi-fluid stability problem. For a linear interfacial perturbations, the commonly used method is to extend the basic flow of one fluid into the other fluid domain in considering perturbation equations. However, when the instability is nonlinear this seems inappropriate to do. 4.2 Three-dimensional nonlinear stability bounds We turn now to the analysis of (4.9), which we can write as: [TW(U + U)-TI?(U)}. (Term A) (Term B) (Term C) (Term D) (4.10) and examine each of terms A-D separately. 51 Chapter 4. The Reynolds-Orr equation and Stability bounds 4.2.1 The Bingham fluid region, n2(t) We treat first the Bingham fluid region. This is an annular region surrounding the Newtonian fluid region, fii(i), which is a cylinder of radius R4, centered at (xc(t),yc(t)). By our earlier assumptions, fii(£) is surrounded by a region of Bingham fluid of finite width, within which both the basic flow and the perturbed flow are unyielded. To formalise this, we assume that there exists h > 0 such that: 7 i » = iij(U) = 0 : V(x ,y) : R* < [(x - xcf + (y - y c ) 2 ] 1 / 2 <Ri + h. (4.11) We denote by ffe.pCO, this ring of unyielded plug surrounding fii(i), see Figure 4.1. Note that there may be other unyielded fluid outside this ring, in f^lf^.p- Since we consider an infinitely long pipe, rotational movements of the plug region are only possible about axes that are parallel to the z-axis. The plug is treated as a rigid body. Thus, at time t, the plug motion can be described as combination of a linear motion, say (uc,vc,wc + Wp), and a rigid body rotation about an axis through (xc(t),yc(t)), parallel to the z-axis. Recall that (uc,vc) is the speed of {xc(t),yc(t)), and here Wp denotes the speed of the plug region in the basic flow. Let u*: u* = {u*,v*,w*) = (uc, vc, wc) + (~(y- yc)uc, (x - xc)ibc, 0), (4.12) where uc(t) denotes the instantaneous rotation of the plug about the z-axis through (xc't),yc{t)). Thus, for any x G Cl2,p(t), u* + U, gives the perturbed velocity of the fluid in the rigid plug, i.e. u* = u. Following the analysis in [10], by considering pointwise the integrand in term A of (4.10), we have: (TermA) = -j^ J f r b » [ 2 | ( C + t i ) - ryI2l(t/)] dx, £ -^Lj^dx-i^b {uUx- (413) This last term, on using the divergence theorem becomes: 52 Chapter 4. The Reynolds-Orr equation and Stability bounds Figure 4.1: Geometry of the perturbed flow. The boundary integrals vanish on the wall, where U{ = 0, and on the ends of the region, due to periodicity in z. At the interface, we use (4.12): I dUi [2] m-K—m /L/2 r2n / [(uc — RiCjc sin 8) sin 9 — (vc + RiUJc cos 8) cos 8] d8 •L/2 Jo = 2nL(RiUc)2. (4.15) The dissipative term A in (4.10) therefore satisfies the following bound: 27rL ( i? i w c ) 2 . (Term A) < ^ / j2 dx = - —^[ [ | | ^ - | 2 da; Rem Thus, rotation of the plug region increases the dissipation, as is to be expected. (4.16) We consider now the inertial term (B) in (4.10). We note first that — 75y = 0 if x <E f22,P Assuming that bounds can be defined of form: max max (x,y)eii; \(xc,yc)\<a = \vcwi (4.17) (4.18) These constants will depend on the parameters (x^y1)?, B, m, Ri. 53 Chapter 4. The Reynolds-Orr equation and Stability bounds To bound uc and vc, we use (4.12) to write: 2uc = u(x,yc-((Ri + h)2-(x-xc)2)^,z) + u(x,yc + ((Ri + h)2-(x-xc)2)^,z), -(Ri + h) <x-xc<(Ri + h), (4.19) 2vc = v(xc - {{Ri + h)2 - (y - yc)2)^,y, z) + v(xc + {{Rt + hf - (y - yc)2)^,y, z), -(Ri + h) <y-yc< (Ri + h), (4.20) i.e. we use the values of u, v on boundary of the unyielded plug interface, to cancel out rotational terms. Writing, u(x, yc - ((Ri + h)2 - (x - xc)2) 2, z) = / J - l - x -yc-((Ri+h)2-(x-xc)2)^ gu — (x,y,z) dy, 2 dy and an analogous expression for u(x,yc + ((Ri + h)2 - (x - xc)2)?,z), we integrate (4.19) with respect to x, to give the following bounds: ,du, , I — I / du2 2uc(Ri + h)L < f \^\dx,< [ & d x , Jn2r dy Jn.ino.r, dy < [n(l-(Ri + h)2)L]^2 j dx 1/2 (4.21) 'n2|o.2,p dy where 02,c is as illustrated in Fig. 4.2, and we have used the Cauchy-Schwarz inequality. Simi-larly, we have: r i V2 / l - l Jn2\n2„ dx 2vc(Ri + h)L < [it(l-(Ri + h)2)L]ll2 dx (4.22) Now let AB,I, AB ,2, be defined as follows: AB,I = sup / n 2 | n 2 p ^ + s ) * d a ; s u p ~7 \dv±\2 A—: uevB,0 io.2\Q.2J dxA a x AB,2 = SUp 'n2|n2,p' dxj In2 ™2 d x uevBfiIU2\^\2 dx (4.23) (4.24) The test space VB,Q above is the closure of VB,O> which consists of functions that satisfy the following conditions: 54 Chapter 4. The Reynolds-Orr equation and Stability bounds Figure 4.2: The domain fi2,c-1. U G C ° ° ( Q 2 ) , 2. it, is solenoidal, i.e. | | * = 0, 3. 7y(ii) = 0 for x £ £1%P, 4. u = 0 on the walls of the pipe: x2 + y2 = 1. 5. u is periodic in 2, with period L . We consider Vg,o to be the closure of VB,O with respect to the [/^(f^)] 3 norm. Combining all the above bounds, terms A and B of (4.10) satisfy: 2nL(RiL)c)2 (Term A ) + (Term B) < 1 1 f ,dui_ 'dx HePl (4.25) where IV7U/I A • ^W\max(l - (Ri + fr)2)1^! - g j ^ A g = | | V ^ | | W A B , I + ^ ^ - ^ . (4.26) In §4.4 we discuss (4.23) and (4.24) further. 55 Chapter 4. The Reynolds-Orr equation and Stability bounds 4.2.2 The Newtonian fluid region, fii(t) We turn now to terms C & D in (4.10), dealing first with the dissipative term. Consider the velocity field u', defined by: u' = u-u*. (4.27) We note the following: 1. Since u* is the perturbation of the plug velocity in Q,2,p and the velocity is continuous across the interface, at [(x - x c ) 2 + (y - yc)2]1^2 = Ri, it follows that: u' = 0, xedflxit). (4.28) 2. Since u* consists of only a linear translation and a rigid body rotation: 7fc(u*)=0, => 7 f c ( u ) = 7 « ( u ' ) - t 4 ' 2 9 ) Using the above, the dissipative term in (4.10) on fii(i) is: (4.30) We now introduce (4.27) into the inertial term in (4.10) on fii(t): f dW dW dW dW (Term D) = - J^(u' + u*) — + W + + —+ ^{w + wc) dx f . ,aw ,dW. , r ,dW bWs ,dW aw„ , , f r , x . dW . .„ dW, , . ,dW ,dW, , - [(y - VcPc-te - - xo)^]w +in^+v -^]we dx r . ,dw aws ,aw dw., , - / [uc(-£- + + « c ( " a - + «—)H da; ax dx c 9y aj/c /" dW dW ~ J [(y~ vJ"c~far ~(x~ ^ c ) w c — ] u i c da; (4.31) In order to bound the above expressions, we need first to find bounds for both wc and Cbc. At time t we define a cylindrical coordinate system (r',61) , centred at (xc(t),yc(t)), see Figure 56 Chapter 4. The Reynolds-Orr equation and Stability bounds Figure 4.3: Coordinate system (r',6'). 4.3. Now wc is the perturbation of the axial component of the plug velocity, and is constant. In terms of the {r',6') coordinate system, the outer wall of the pipe is denoted by r' = r0(8'). We have the following: cL/2 r2w / wc(8') dd'dz •L/2 JO rL/2 m fRi+hQ = (Ri + h) / / g d b W d * J-L/2 JO Jr0(6') Pr' fL/2 /°2TT «•„(«') Q < (Ri + h) / / / J-L/2 JO JR dr'de'dz I Ri+h run rlix proW) Qw f < / / / \^j\r'dr'd8'dz = J-L/2 Jo JRi+h dr J n f!2 1^ 2,p dec < [ T r a - ^ + ^ j L ] 1 / 2 / \ — I 1,1/2 dec (4.32) 57 Chapter 4. The Reynolds-Orr equation and Stability bounds To bound uc we again use (4.12) to write: 20c((Ri + h)2- (x-xc)2)* = u(x,yc-{(Ri + h)2-(x-xc)2)^z) (4.33) -u(x,yc + ({Ri + h)2 - ( x - z c) 2)5, z), -(Ri + h) < x - xc < (Ri + h), 2uc((Ri + h)2 - (y - yc)2)^ = v(xc + ((Ri + h)2 - (y - yc)2)Ky, z) (4.34) -v(xc - ((Ri + h)2 - (y - yc)2)Ky, z), -(Ri + h) <y-yc< (Ri + h). As before, we write u(x,yc - ((Ri + h)2 - (x-xc)2)^, z) and u(x,yc + ((Ri + h)2 - (x - xc)2)^, z) in terms of integrals of and integrate (4.33) between xc ± (Ri + h). We perform similar operations on (4.34), u>cirL(Ri + h)2 < f Adx< [ Adx, (4.35) Jn2,c dy Jn2\n2,„ dy ucnL(Ri + h)2 < / | ^ | dx. (4.36) Jn2\n2.v dx 'n 2|n 2, P Finally, combining the above expressions and using the Cauchy-Schwarz inequality: 1/2 ( l - ( R i + h)2fl2 lUJcl ~ (27rLy/2(Ri + h)2 j | ^ i | 2 dx n2\n2,P dxj (4.37) We now consider upper bounds for the individual terms in (4.31). We first note that if fti is surrounded by a region of unyielded Bingham fluid, then the basic flow W is such that W = W(r'), and furthermore, for x G Q\: W(r') = WN,B l - ( ^ r ) 2 Hi + W r ^ To see this, consider the problem for W(x,y;xc,yc) in (see §3.1). If the Bingham fluid is unyielded then W = Wp = constant. The function W - Wp satisfies Poisson's equation with homogeneous Dirichlet boundary conditions on a circle of radius Ri, which has a unique solution W — Wp = h(r'), i.e. is axisymmetric. 58 Chapter 4. The Reynolds-Orr equation and Stability bounds Now let AAT,I, Ajyt2) Ajv,3, be denned as follows: /n t r'ur'W dx JQl[r'w}2 dx ANtl = sup 7 1 , (4.38) uevN,0 Jni & d x A",4 - .sup r U ; ^ - (4-41) A " ' 2 s - S U P f Ska H T ' ( 4 ' 3 9 ) fn w2 dx AiV'3 S /Vp da;' (4'40) JUl[r'ur'}2 dx where ur> is the radial component of it, i.e. in the direction of r'. The test space VN,O above is the closure of VJV,O, which consists of functions that satisfy the following conditions: 1. u € C°°(Oi), 2. u is solenoidal, i.e. = 0, 3. u = 0 at the boundary of the Newtonian region: (x — xc)2 + (y — y c ) 2 = Tj. 4. i i is periodic in z, with period L . We consider VN,O to be the closure of V/v,o with respect to the [ i / 1 ^ ) ] 3 norm. In §4.4 we discuss (4.38)-(4.41) further. 59 Chapter 4. The Reynolds-Orr equation and Stability bounds Finally, we have for (4.31): (TermD) < f | H 2 dx + + + [2A^, 2 (1 - (Rj + / i ) 2 ) ] 1 / 2 Ri(Ri + h) ^ | V c ^ | m a x ^ [ 2 A w , 3 ( l - (Ri + fr)2)]1/2 Ri + h WN .max [kNA(l-(Ri + h)2)}ll2 Ri(Ri + h) .max (I-(Ri + h)2) \ l l l ^ l 2 dx 1/2 " h 9XJ l l ^ l 2 dx" 1/2 ' .7r /" i duj I n 1/2 dx n 2 | n 2 , p 1/2 dx l l ^ l 2 dx 1/2 " 7r 7n2|n2,„ n 1/2 dx / JO, 3\/2(^ + /i)2 7r / 2 ,? |V c W| m a x( l - ( i i i + ^ ) /" 7f! 0 U j n 2 | n 2 , p i ,2 dx 2 02|O.2,p 9x3 We now 2 v / 2(i^ + h)2 combine the above bounds for terms C and D in (4.10): (Term C) + (Term D) < AN j | | ^ | 2 dx + CN j | | ^ | 2 7n 2 |n 2, p dxj JUl dxj " / (4.42) | duj l2 Jn2\n2,P dxj +2BN m dx l l ^ l 2 dx 1/2 " .A 'Rem f 1 ^ 7« i J 1^ 2)^ 2, p | ^ | 2 dx 1/2 I2 dx (4.43) (4.44) where Ajv, -Bjv and CAT are the positive parameters: irR2(l-(Ri + h)2) r -Bjv = Civ = V2(Ri + h)2 3Ri 2 2WAT ; m a x A j y ; i R2 (4.45) (4.46) (4.47) 60 Chapter 4. The Reynolds-Orr equation and Stability bounds 4 . 3 The stability bound We combine (4.10) with the bounds (4.25) & (4.44)to give 2DiiU d X * { A N + B N + AB-R^)jJdX-^ d X ReW m , f , du'i +(BN + CN- )( \ p - ? d x . (4.48) Suppose the following 2 conditions hold, (for nonlinear stability): *•* < AM + Bpj + AB' ' <4'49) m i j j v + CN Assuming (4.49) & (4.50) hold, we rewrite (4.48) as: ,2 l^i f ^ d x ^ -(-^ ~ AN - BN - AB) f j(u)2 dx - (AN + BN + Ab)2TTL(R4QC)2 2 Dt JN 2 ReW JQ2 ~{RW\ ~ Bn~ Cn) L ^ dX' (4' 51) < - A f 7 (tt) 2 dec, (4.52) < - cA f u2 dx, (4.53) Jn where 1 m A = ram{(-^ -AN-BN- AB), (^ - BN - CN)} (4.54) In moving between (4.52) & (4.53), we have used Korn's inequality on Q, with c the result-ing constant. Hence, the following energy stability bound holds for the exponential decay of |2 l u l iL 2 (n)-\u\\h{n)(t) <\\u\\2L2m(0)eW'2cXt. (4.55) Note that, through (4.12), the speed at which |(x c,y c)| changes is related to the plug velocity, which we have shown decays exponentially with t, provided (4.49) & (4.50) hold. We can thus use the decay bound (4.55) to show that |(x c,j/ c)| remains bounded for all time. 61 Chapter 4. The Reynolds-Orr equation and Stability bounds 4.3.1 The bound for \(xc, yc)\ We have assumed throughout that |(a;c,yc)| < a, which we now verify and also determine an estimate for a. First note that (xc, yc) satisfies the following initial value problem: jt(xc,yc) = (uc,vc), (xc,yc)(0) = (0,0). (4.56) To bound uc and vc, we use (4.12) again, and note that for r' € (Ri, Ri + h): 2uc = u(x,yc-((r')2-(x-xc)2)5,z) + u(x,yc + ((r'f-(x-xc)2)Kz), -r' < x - xc < r', (4.57) 2vc = v(xc-((r')2-(y-yc)2)Ky,z) + v(xc + ((r')2-(y-yc)2)^,y,z), -r' <y-yc< r', (4.58) we now integrate over 0,2,p'-1/2 udx < (Ln[(Ri + h)2 - r f ] ) 1 / 2 | / u2dx\ and similarly for vc, Then using (4.55) cLTr[(Ri + h)2 - R2} = [ udx< (LTr[(Ri + h)2 - r2])1'2 f u2dx i \n\ \\U\\L2(Q)(0) _cXt * (L,{(Ri + h ) 2 - R 2 } y f 2 ^ ' ( 4 - 5 9 ) i ^ H u l lL2 (n ) (0) _ c X t l V c m s (L.[(Ri + h ) 2 - R 2 ] y / 2 ^ ( 4 - 6 0 ) We now integrate (4.56) with respect to t and insert the above bounds on |uc|(t) and |uc|(i), to give N(*)<yo M ( ^ < c A ( M ( ^ + , )2_ i,2])1/2> ^^^L^AY-R2])^-(4.61) Thus, finally we have: l ( X c ' y c ) l ( i ) * cX(L,[(Ri + h)2-R2]y/2 = ^ ( 4 6 2 ) 62 Chapter 4. The Reynolds-Orr equation and Stability bounds 4.4 Parametric variations in the stability bounds The stability bounds on the Reynolds numbers (4.49) & (4.50) which lead to the decay bounds (4.55) & (4.62) have been arrived at after considerable analysis and contain a number of constant parameters defined en-route in §4.2.1 and §4.2.2. The first point to establish is whether or not these parameters are well-defined. Secondly, if possible, we would like to give some approximate numerical value of the constants. Thirdly, we would like to establish how the bounds and the decay constants vary with each of the process parameters that define the basic flow, (xc,yc), B, m, and Ri. 4.4.1 \VW\2,max, | V c W | m a a ; , Wjv,maa;. Since we consider only bounded \(xc,yc)\ < a, for a given by (4.62), and our basic solutions depend on (xc,yc), B, m, and Rt, the constants |VW|2,moa:, | V c W | m o x , WN,max are well defined provided that W(x, y\xc, yc) is differentiable with respect to each of the above parameters. We have made no effort to determined values for any of these constants. However, it is clear that such values could be approximated by using the numerical solutions of §3.3. Similarly, parametric dependencies could be explored by using the numerical solutions. From inspection of our numerical solutions, we believe that Wjv,max increases with both B and Ri, but decreases with m. The parameter |VVF|2,max is generally 0(1), it increases with m, decreases with Ri and may either increase (small B) or decrease (large B) with B. The latter effect is explained physically, since for small B the yielded region may decrease with little variation in the plug velocity, whereas larger B eventually results in the Bingham fluid becoming stationary, when | V W | 2 , ,max — 0- To compute |V cW^m ax is numerically a time consuming problem. We note that W = 0(1) since this is the mean value, and that since the flow rate constraint is always satisfied, changing (xc,yc) will simply increase W in some parts of the domain whilst decreasing W in other parts. It is therefore thought that IVcWImax is numerically of 0(1) and fairly insensitive to changes in B, m, and Ri. 63 Chapter 4. The Reynolds-Orr equation and Stability bounds 4.4.2 A B ) 1 , A B ) 2 . Constants such as A ^ , i and AB,2 can be computed from an eigenvalue problem, via the Euler-Lagrange equations for the appropriate minimisation problems. The constant parameters are found as the smallest positive eigenvalues. In the case that (xc,yc) = (0,0), for all t, (i.e. the Newtonian region does not move), the constant A#,2 does not appear and the constant AB, I has a slightly simpler form. This is discussed briefly in §5.1.1. In general, it is very hard to progress with sharp bounds for (4.23) and (4.24) when the asymmetric case is considered, due to the shape of Cl2\Cl2,p. However, it is possible to establish that A#,i and AB,2 are well defined, as follows. We know that, \w < \ JRi+h dw dr' dr' < r / r'5 i r (Ri + h) 2 JRi+h dw dr' dr'. Now, f w2dn2\tt2iP < UL r'w(e'] r'w(9') I Ri+h (Ri + h) Using Cauchy-Schwarz inequality we get JR., +h dw dr" dr" dr'dO'dz f w*dn2\n2,p < f [l^M^[r'w(e>)-(Ri + h)}2 Jn2\n2iP J J (iti + n) JRi+h Applying the Mean value Theorem for Integrals we have, / w " d x * 2^R:+h,h){r'w(e'*)-(Ri + h)? [ n2|n 2,p but Jn2\n2,p \ dxj J X Jn2\n2„ V dr" ) dx dw dr" dr"dd'dz dw dr" dx (4.63) Therefore, we have that: 64 Chapter 4. The Reynolds-Orr equation and Stability bounds In order to bound KB i we use the result above and the Cauchy-Schwarz inequality, It / f i 2 | n 2 , P ( m$ ) d x \fn2\n2,P{^) dx) \(-/n 2 |n 2 , p ) d x uwdx f iv^dx J f 2 2 | n 2 i P Si n 2 | n 2 , uwdx \ \/fi2|fi2,p ( | ^ | ) Therefore, finally we have that: < i n 2 | o . 2 , p ( ) d a ; This establishes that A^, i and KB,2 are well-defined. Its quite clear that our bounds are very conservative. The chief approximations used are the following: (i) Cauchy-Schwarz, (ii) Effectively including all velocity components in our bound, (iii) Complete neglect of the divergence-free constraint, (iv) Complete neglect of the boundary conditions at the plug. To investigate the parametric dependency of KB,i and KB,2, we examine our numerical results in §3.3. The factor [r'w(0'*) - (Ri + h)]2 which appears above is effectively related to the size of the unyielded region, as is R4 + h. We tentatively conclude that Ks,k will increase with Ri and will decrease with B. The effects of changing m are unclear. 4.4.3 AJV,I , Ajy,25 AJV,3? Ajv,4-Estimates for AJV,I, Ajy,2, KN,3, AJV,4 are much easier due to the circular domain and Dirichlet conditions. Each bound can be written as a eigenvalue problem for a different set of Euler-Lagrange equations. These eigenvalue problems are solved by first substituting a normal mode and then computationally solving the resulting parametric eigenvalue problem. In principal we could do this to get numerical values. This is however quite time consuming to do, and especially 65 Chapter 4. The Reynolds-Orr equation and Stability bounds since the minimal eigenvalues are not always found to lie at the boundaries of the parameter regime, (i.e. meaning that the 2D axisymmetric or azimuthal problems often do not yield the minima, and hence a 3D parametric problem must be solved and explored numerically). Instead we just establish that these bounds exist and then give an upper estimate for the bounds. For simplicity we map the Newtonian region into the unit circle, which simplifies the dependency on Ri. The bound for A^, i is due to Joseph [7]: AJV,I< ^ 81.49' This bound is calculated from numerical solution of the following eigenvalue problem: iarDLf - (N2 + a2r2)Lw - 2a2rDw = pN2f p[(N2+ a2r2)w + iarDf] = L2f + 4iaLu L = -D(rD) - - j -/ = Df = w = 0 at r — 1 f,w bounded at r — 0 where f = z,N = l,a = 1.07 and p = -fa. For Ajy^) Ajv,3, AJV^, we relax the constraints in the definition of these bounds, as follows: /o \f'w}2 dx A n> 2 - - s u p r<4 dx' ( 4 - 6 4 ) f vP dx Ajv'3 - -T r f e - (4-65) f0 \r'urA2 dx A n > 4 ~ - s u p r Uv2, (4'66) Note that by eliminating the additional terms in the denominators, we are effectively throwing away the divergence-free condition. After writing (4.64) as an eigenvalue problem and solving it, we find that the least positive root of the Bessel function Jo(^) gives the desired value. Note that p = and thus for A/v,2 we have: AJV,2 < R * 23.13275 66 Chapter 4. The Reynolds-Orr equation and Stability bounds Following the same procedure as above with (4.65), p = is found the least positive root of Jo(y/p). The bound for Ajv3 is A j V ' 3 - 5 J 8 3 1 9 ' We can see that (4.66) is equivalent to (4.64), thus A w ' 4 - 23.13275-Again we note that our bounds for Aj\r,2i AJV,3, Aj\r,4 are likely to be conservative. We expect that the numerical values in the denominators above will be of O(102), (by analogy with Joseph's results). 4.4.4 Variation of (4.49) & (4.50) If we now try to translate the parametric understanding above into the definitions of AN, BN, CN and As, which effectively define the stability bounds (4.49) & (4.50), we are unfortunately defeated in making any statements of any value! Changing any of the process parameters B, m or Ri does not definitively move all of these parameters in one direction or the other, so we can not definitively make any statement regarding the effects on stability. 67 Chapter 5 Discussion and conclusions 5.1 The bounds (4.49) & (4.50) This thesis is concerned primarily with deriving (4.49) & (4.50), which has been quite a hard problem. Due to the complexity of the derivation and the many functional analytic approxi-mations made, our eventual bounds really just say that the flow is stable for sufficiently small Reynolds number. A similar statement can be made for an arbitrary single phase flow, but not for a multi-phase flow. It is important to realise that the structure of our (multi-phase) flow is preserved under our stability bounds. This is a great step forward, since it implies that we may actually be able to observe this flow experimentally. Indeed, given the general nature of our flow, (i.e. 2 fluids which may mix, finite periodic perturbations subject only to bounds on the shear stress perturbation), the type of stability bound we have is as good as we might expect. In order to make the motion of the Newtonian region as small as is possible, by inspecting (4.62), it is possible to do one of 3 things: (i) Increase the size of the plug region (ii) Increase the decay constant A, by moving away from (4.49) & (4.50) as far as possible, (iii) Limit the initial velocity perturbation. If we experiment with (4.49) & (4.50), using our estimates from §4.4, we usually find that stability is achieved at Reynolds numbers of 0(1). This is not very dramatic. However, actual stability is likely to be much better than predicted by (4.49) & (4.50). Our approach throughout has been to make (sometimes quite crude) approximations and in any case the energy method 68 Chapter 5. Discussion and conclusions is not an exact method. Much of the approximation could be improved, but in the end, to do this requires intensive computation, (e.g. finding exact values for A ^ i and A ^ ) . Whether this is worthwhile to do for a general problem is debatable. In this context, the important thing is to have achieved is a nonlinear bound that preserves the structure of the flow. 5.1.1 Stability where \(xc,yc)\ = 0 For practical purposes, we might try to get a better estimate of sharper stability bounds by considering a simpler problem, namely that when \(xc,yc)\ = 0. Movement of the plug away from its central position only occurs when there is a net force on the yield surface in one particular direction. In a well-controlled process, one would try to limit asymmetry in the stress field and velocity field, as well as bounding the magnitude of any unwanted perturbations from the axisymmetric base case. This might be quite hard to do practically, but supposing it could be done, the best case that we are likely to encounter is when \ {xc, yc)\ = 0, which we can briefly consider here. For \(xc,yc)\ — 0, much of our analysis is greatly simplified. The Reynolds-Orr equation is now: We can proceed as before, but now we will only need to define and use AJV,I and in place of AB,I we will have to compute only f n i n W'(r)urw dx A B = sup — . (5.2) uevBi0 Ju2\u2Js^r dx The Bingham domain is now a concentric annulus and the above problem is at least tractable. By following through the analysis it becomes clear that our Reynolds number bounds will be significantly improved, but we have not computed any numerical values. We have made some efforts to numerically solve the eigenvalue problem associated with (5.2). Note that this eigenvalue problem is over specified in that it has too many boundary conditions at the yield surface, via the definition of VB,O- This problem is also quite similar to that studied in [10], where only approximations were made. It is not clear what effect having too many boundary conditions will have, and our numerical solutions have not yet resolved this problem. 69 Chapter 5. Discussion and conclusions 5.1.2 One-dimensional stability Another area in which some test computations can be carried out is in considering one-dimensional stability problems. Either a purely axial or azimuthal perturbation could be con-sidered. We can see from the Reynolds-Orr equation that such perturbations, (which necessarily impose |(x c ,y c)| = 0), are nonlinearly stable. The point of doing this type of computation is to get an idea of what would be the best possible decay rate for a stable perturbation. 5.2 Summary of the main contributions of the thesis We have considered an axisymmetric basic flow in which a Bingham fluid surrounds a Newtonian fluid. The flow rates and rheologies are such that the Bingham fluid is unyielded at the interface. The thesis has considered the nonlinear stability of this flow. The contributions are as follows. • We demonstrate that by using a visco-plastic fluid, (i.e. a fluid with a yield stress), as the lubricating fluid in the above multi-layer pipe flow, we eliminate interfacial instabilities, i.e. by preserving an unyielded region at the interface which cannot deform. • We derive a conditional nonlinear stability bound for a quite arbitrary 3-dimensional perturbation. Under these bounds the structure of the flow is preserved. • Our bound allows motion of the plug and Newtonian region, but we show that this motion can be bounded for all time. • We compute 2-fluid basic flows using an augmented Lagrangian method, that accurately represents the unyielded regions of the flow. 5.3 Future work There are many possible avenues for further work, and we list a few. • In demonstrating weak nonlinear stability of these flows, we motivate the possibility of demonstrating these flows experimentally. This is one avenue for future research. Once demonstrated experimentally, this would open the way for a range of applications to 70 Chapter 5. Discussion and conclusions industrial processes where multi-layer processing might be of interest, e.g. coated wires and tapes, co-extruded laminates, etc. • It would be interesting to compute these transient flows directly, in 2D or 3D, probably using an augmented Lagrangian approach. • The eigenvalue problem (5.2) is related to many others that occur in stability problems involving Bingham fluids and it is interesting in that it is over-specified. A deeper study of such problems could be interesting. • An exploration of the simplified problems outlined in §5.1.1 and §5.1.2 would be worth-while. 71 Bibliography [1] M . Fortin and R. Glowinski. Augmented Lagrangian Methods. North-Holland, 1983. [2] L A . Frigaard. Super-stable parallel flows of multiple viscoplastic fluid. J. Non-Newt. Fluid Mech., 100, pages 49-76, 2001. [3] L A . Frigaard, O. Scherzer. Uniaxial exchange flows of two Bingham fluids in a cylindrical duct. IMA J. Appl. Math., 61, pages 237-266, 1998. [4] C . E . Hickox. Instability due to Viscosity and Density Stratification in Axisymmetric Pipe Flow. Phys. Fluids, 14(2), pages 251-262 1971. [5] E . J . Hinch. A note on the mechanism of the instability at the interface between two shearing fluids. J. Fluid Mech., 114, pages 463-465, 1984. [6] A.P. Hooper and W . G . Boyd. Shear flow instability at the interface between two viscous fluids. J. Fluid Mech., 179, pages 201-225, 1987. [7] D.D. Joseph and S. Carmi. Stability of Poiseuille Flow in Pipes, Annuli, and Channels. Quart. Appl. Math., 26, pages 575-599, 1969. [8] D.D. Joseph, Y . Y . Renardy and M . Renardy. Instability of the flow of two immiscible liquids with different viscosities in a pipe. , J. Fluid Mech., 141, pages 309-317, 1984. [9] D.D. Joseph and Y . Y . Renardy. Fundamentals of Two-Fluid Dynamics Vol. I and Vol. II. Interdisciplinary Applied Mathematics, Springer, 1993. [10] C. Nouar and L A . Frigaard. Nonlinear stability of Poiseuille flow of a Bingham fluid: theoretical results and comparison with phenomenological criteria. J. Non-Newt. Fluid Mech., 100, pages 127-149, 2001. [11] A. Pinarbasi and A. Liakopoulos. Stability of two-layer Poiseuille flow of Carreau-Yasuda and Bingham-like fluids. J. Non-Newt. Fluid Mech., 57, pages 227-241, 1995. [12] S.G. Yiantsos and B.G. Higgins. Linear stability of plane Poiseuille flow of two superposed fluids. Phys. Fluids, 5i,pages 3225-3238 1988. [13] C.-S. Yih. Instability due to viscosity stratification. J. Fluid Mech., 27(2), pages 337-352, 1967. 72
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Nonlinearly stable multilayer viscoplastic flows
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Nonlinearly stable multilayer viscoplastic flows Moyers Gonzalez, Miguel Angel 2002
pdf
Page Metadata
Item Metadata
Title | Nonlinearly stable multilayer viscoplastic flows |
Creator |
Moyers Gonzalez, Miguel Angel |
Date Issued | 2002 |
Description | When two purely viscous fluids flow in a parallel multi-layer shear flow down a pipe or other duct, the flow is frequently unstable. The primary instability arises at the interface, which is linearly unstable to certain wavelengths. Physically, for any perturbation of a planar interface, the less viscous fluid accelerates the more viscous fluid; see [5]. Although it is possible to suppress short wavelength instabilities through surface tension and also to achieve a more stable flow through arranging the viscosities of the fluid layers, (see e.g. [9]), none of these methods actually eliminate the underlying interfacial instability. In this thesis we demonstrate that by using a visco-plastic fluid, (i.e. a fluid with a yield stress), as the lubricating fluid in a multi-layer pipe flow, we do actually eliminate interfacial instabilities. This fact has been established recently in [2] using a linear stability analysis. In this thesis the analysis is extended to nonlinear stability, applying the methods developed in [10]. Our basic flow is an axisymmetric flow in which a Bingham fluid surrounds a Newtonian fluid. The flow rates and rheologies are such that the Bingham fluid is unyielded at the interface. We consider the nonlinear stability of this flow. This nonlinear stability problem is interesting for a number of reasons. Firstly, for a nonlinear perturbation it is possible for the unyielded plug region, (and the Newtonian region that it surrounds), to be perturbed and move. For the linear stability problem in [2], such motions do not occur since the stress perturbations on the interface and yield surface are linear and periodic, (via a normal mode approach). Thus, they integrate to zero over the volume of the plug and there is no motion. To allow for motion of the plug and Newtonian fluid region, whilst preserving a stable interface, is the main challenge of this thesis. The stable interface is preserved by bounding the perturbation in the deviatoric stress, which ensures that an unyielded region will still surrounds the Newtonian region. To carry out the analysis for a moving region we then consider a perturbation about a basic flow that is asymmetric. We find stability bounds, (conditional on the stress perturbation), below which the velocity perturbation from the asymmetric basic flow decays. We then show that the displacement of the plug region is also bounded for all time. Thus, in the above way we prove that the axisymmetric basic flow is stable but not asymptotically stable, i.e. conditional nonlinear perturbations will decay to an asymmetric solution that is close to the axisymmetric basic flow. Secondly, the asymmetric basic flows are themselves interesting. These flows fall into a class of flows that have been studied in [3], from the point of view of proving existence and uniqueness, in the context of examining various problems in cementing. Here we compute these flows using the augmented Lagrangian approach. This method is special in that it manages to compute the unyielded behaviour exactly in regions where the stress is below the yield stress. This is the first time, to our knowledge, that this method has been used to compute multi-fluid flows of visco-plastic fluids. Thirdly, the flows we study are special in that the structure is preserved, (i.e. due to one of the fluids having a yield stress). For any multi-phase fluid problem this is unusual. It is really this preservation of structure that has allowed the application of nonlinear stability methods. As such, this study is one of very few multi-phase flow problems that we know of, where an energy method has been effectively applied. Finally, in demonstrating weak nonlinear stability of these flows, we open the way for a range of applications to industrial processes where multi-layer processing might be of interest, e.g. coated wires and tapes, co-extruded laminates, etc. |
Extent | 3215873 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-09-28 |
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.0080079 |
URI | http://hdl.handle.net/2429/13290 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2002-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2002-0505.pdf [ 3.07MB ]
- Metadata
- JSON: 831-1.0080079.json
- JSON-LD: 831-1.0080079-ld.json
- RDF/XML (Pretty): 831-1.0080079-rdf.xml
- RDF/JSON: 831-1.0080079-rdf.json
- Turtle: 831-1.0080079-turtle.txt
- N-Triples: 831-1.0080079-rdf-ntriples.txt
- Original Record: 831-1.0080079-source.json
- Full Text
- 831-1.0080079-fulltext.txt
- Citation
- 831-1.0080079.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0080079/manifest