B U B B L E P R O P A G A T I O N T H R O U G H VISCOPLASTIC FLUIDS by NEVILLE DUBASH BMath (Applied Mathematics) Hons. Co-op, University of Waterloo, 2001 A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF THE REQUIREMENTS FOR THE D E G R E E OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Mechanical Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA • August 2003 © Neville Dubash, 2003 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 reference 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 Mechanical Engineering The University of British Columbia Vancouver, Canada Date II A<-u^u.s Abstract In this thesis we consider the propagation of an air bubble in a cylindrical column filled with a viscoplastic fluid. Because of the yield stress of the fluid, it is possible that a bubble will remain trapped in the fluid indefinitely. We restrict our focus to the case of slow moving or near-stopped bubbles. Using the Herschel-Bulkley constitutive equation to model our viscoplastic fluid, we develop a general variational inequality for our problem. This inequality leads to a stress minimization principle for the solution velocity field. We are also able to prove a stress maximization principle for the solution stress field. Using these two principles we develop three stopping conditions. For a given bubble we can calculate, from our stopping conditions, a critical Bingham number above which the bubble will not move. The first stopping condition is applicable to arbitrary axisymmetric bubbles. It is strongly dependent on the bubble length as well as the general shape of the bubble. The second stopping condition allows us to use existing solutions of simpler problems to calculate additional stopping conditions. We illustrate this second stopping condition using the example of a spherical bubble. The third stopping condition applies to long cylindrical bubbles and is dependent on the radius of the bubble. In addition to our stopping conditions, we determine how the physical parameters of the problem affect the rise velocity of the bubble. We also conduct a set of experiments using a series of six different Carbopol solutions. From the experiments we examine the dependence of the bubble propagation velocity on the fluid parameters and compare this to our analytic results. We find that there is an interesting discrepancy for low modified Reynolds number flows wherein the bubble velocity increases with a decrease in the modified Reynolds numbers. We also compare our three stopping conditions with the data. It appears that all the stopping conditions seem to be valid for the range of bubbles examined despite the fact that when applying the second and third stopping conditions most bubble shapes are not well approximated by a sphere or a cylinder. u Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Acknowledgement Chapter 1. 1.1 1.2 1.3 3.4 4.2 4.3 Problem Description 8 8 11 Variational Formulation 15 Stopping Condition Results 15 15 18 18 20 23 25 25 31 First Stopping Condition 4.1.1 Condition On B For No Motion 4.1.2 Surface Integral Term 4.1.3 Bubbles Which Never Move 4.1.4 B o u n d i n g - j f ^ T Second Stopping Condition 4.2.1 Second Stopping Condition for a Spherical Bubble Third Stopping Condition: for Long Cylindrical Bubbles Chapter 5. 5.1 5.2 5.3 1 3 6 Introduction Simplification of the Full Problem Variational Inequality (Rate of Strain Minimization) 3.3.1 A Preliminary Inequality 3.3.2 Derivation of the Variational Inequality 3.3.3 An Alternate Formulation 3.3.4 Existence/Uniqueness Stress Maximization Principle Chapter 4. 4.1 1 Physical Setup Non-dimensionalized Equations Chapter 3. 3.1 3.2 3.3 Introduction Viscoplastic Fluids Previous and Related Work Outline of Thesis Chapter 2. 2.1 2.2 viii 31 31 32 36 3 7 40 41 43 Parameter Dependence 47 Consistency Yield Stress (Bingham Number) Density 48 49 50 iii 5.4 Surface Tension Chapter 6. 6.1 6.2 6.3 51 Experiments 53 Experimental Setup 6.1.1 Experimental Apparatus 6.1.2 Viscoplastic Fluid - Carbopol 6.1.3 Preparation of Carbopol Solutions 6.1.4 Experimental Method Experimental Results Comparison with Analytic Results 6.3.1 Parameter Dependence 6.3.2 Stopping Conditions Bibliography 53 53 54 55 56 57 61 61 62 69 Appendix A . A Result on the Effect of Walls 72 Appendix B . Some Differential Geometry Results for Surfaces 74 Appendix C . Data Extraction and Error Analysis 77 C . l Velocity Calculation C. 2 Shape Dependent Quantities C.2.1 Assumption of Axisymmetry C.2.2 Accuracy of the Edge Detection Method Appendix D . Measurement of Rheological Parameters D. l Yield Stress D.2 Consistency and Power Law Index Appendix E . 77 78 79 79 81 81 81 Optical Distortion Due to Cylindrical Geometry iv 83 List of Tables 6.1 6.2 Fluid properties of the Carbopol mixtures Velocities of bubbles for which •§- > 0.5 in the different Carbopol solutions 57 59 C.l Error in the calculation of shape dependent quantities for test circles 80 v List of Figures 2.1 Physical setup of the problem 4.1 4.2 4.3 4.4 Axisymmetric bubble centred in the column Annular slice of volume, V, over which we are integrating Top view of the horizontal cross section S and Sb Two bubble profiles. For the solid line S = -0.010, and for the dashed line S = 0.012 is nondimensional) Spherical coordinates A long <C L) cylindrical bubble Velocity profile of the flow around the cylindrical section of the bubble 4.5 4.6 4.7 6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 8 Schematic diagram of the experimental setup 54 Curves of velocity versus volume for two Carbopol solutions 58 The complete data set 60 Modified Reynolds number plotted versus Bingham number. For the two solutions the data lie along the same curve (dotted line) 61 Contours of the Froude number, - 7 ^ , plotted as a function of the Bingham number and the modified Reynolds number. The contour spacing is 0.025 62 Froude number as a function of Bingham number for fixed Re* 63 Froude number as a function of modified Reynolds number for fixed B 64 Experimental results plotted with non-dimensional bubble length versus the Bingham number. The solid line represents the curve B = ^^{ + ~ -) d the dotted line represents the curve B = ^ 65 Experimental results plotted with Froude number versus the Bingham number. The dotted line represents B = ^ j . Note that B decreases with increasing Fr, and thus z 6.9 32 33 34 (S 37 42 44 45 z a n faster bubbles are further from the B = ^ line 66 6.10 Experimental results plotted with non-dimensional bubble length versus the Bingham number. The critical Bingham numbers for the data are also plotted (dots). The dashed line is a best fit line through the critical Bingham numbers 67 6.11 Experimental data combined with all three stopping conditions 68 A. l A bubble propagating in two different domains 73 B. l Surface of revolution 75 Cl C. 2 Video frames of an experiment, (a) t = 0 s, (b) t = 1.28 s 77 Profile of a bubble, (a) actual bubble (same as figure C.la), (b) the profile obtained using our edge detection method. The apparent discrepancy in the profile when compared to the video is a result of the cylindrical geometry, for which we have corrected 78 Profiles of a bubble obtained by looking at the left half (dashed line) and at the right half (solid line) 79 C.3 vi D.l D. 2 E. l Strain response after 300 s for a sample of Carbopol at various stresses. The lines are to aid visualization 82 Stress plotted versus rate of strain for a sample of Carbopol. The solid line represents the least-squares fit 82 Schematic of the cylindrical geometry 83 vn Acknowledgement This research was conducted under the supervision of Dr. Ian Frigaard and Dr. Sheldon Green. I would like to thank them both for for their guidance, advice, and encouragement throughout. I would also like to thank the Natural Sciences and Engineering Research Council of Canada, and the British Columbia Advanced Systems Institute for their generous financial support; and the Canada Foundation for Innivation, the British Columbia Knowledge Development, and the University of British Columbia Blusson Fund for their infrastructure funding that enabled the experimental apparatus to be constructed. viii Chapter 1 Introduction In this thesis we examine the propagation of air bubbles through a viscoplastic fluid. Because of the yield stress of the fluid we have the interesting phenomema that bubbles can become trapped in the fluid indefinitely. We are primarily interested in the of case slow moving and stopped bubbles, and conditions that determine whether a given bubble will rise or remain trapped. 1.1 Viscoplastic Fluids Viscoplastic fluids are fluids that can behave as both a viscous liquid and as a rigid solid. Which of these two states the material takes is determined by the stress applied to the material. That is, at low stress the material behaves as a solid, capable of rigid motion and rotation. However, when stresses exceed a certain threshold, called the yield stress, the material behaves as a viscous liquid. The solid and liquid regions, called the unyielded and yielded regions respectively, are separated by a distinct yield surface at which the stress equals the yield stress. Some examples of viscoplastic fluids are mud, cement, lava, and paint. The constitutive models most often used to describe a viscoplastic fluid are the Bingham constitutive model [1] or the slightly more general Herschel-Bulkley constitutive model (2.8)-(2.9). The Bingham constitutive model is a special case (n = 1) of the Herschel-Bulkley constitutive model. In general, foams, slurries, suspensions, and many other industrial fluids behave as viscoplastic fluids. In addition, there are many industrial processes where the generation and behaviour 1 of bubbles play a significant role in the process. In some situations bubbles are undesirable, and while it would be ideal if they could be eliminated, in practice this is not possible. Hence an understanding of how bubbles behave is beneficial in these areas. Some examples of such processes are 1. When drilling for oil or natural gas a situation known as a kick can occur. A kick occurs when the hydrostatic pressure of the drilling mud becomes less than the pressure of the gas in the rock formation. In this case gas from the formation flows into the wellbore and can cause a blowout (where the gas rises up pushing the drilling mud back out the top of the wellbore) [2]. Blowouts should be avoided as they are dangerous to workers, present environmental hazards, and result in a loss of revenue. The drilling mud is generally a viscoplastic liquid. Thus, knowledge of how bubbles propagate could be used to help control kicks and prevent blowouts. 2. A lot of processed foods (e.g., ketchup, mayonnaise, peanut butter) behave as a viscoplastic liquid. In food processing bubbles are unwanted in the processing and especially the packaging steps. Here knowledge of how bubbles are generated and how they propagate could be used to develop methods to prevent or remove bubbles before processing and packaging. Moreover, in other situations bubbles are desired or beneficial to the process. An example of such a process is in plastic moulding. One method of manufacturing bottles is to inject a gas into a mould filled with a molten plastic. This gas forces the molten plastic to take the shape of the mould, and once cooled this becomes the bottle. Again, the molten plastic can be considered to be a viscoplastic fluid. An understanding of bubbles could help to improve and streamline this process. While in this thesis we do not consider any of these problems specifically, the results and analysis, with a little further work, could be applied to these areas. 2 1.2 Previous and Related Work Separately bubbles and viscoplastic fluids have been studied quite extensively and for many years, however, the problem of the propagation of a bubble in a viscoplastic fluid has yet to receive much attention. To the best of our knowledge the problem has only been examined in several specific application and situations. The problems has been examined from the perspective of well drilling, where the gas rise velocity is an important characteristic to know or predict [2] [3] [4]. Santos and Azar [2] measure the gas rise velocities of bubbles rising through an annulus filled with a viscoplastic fluid. They correlate their experimental results with a simple analytic model based on a force balance between viscous drag and buoyancy. To account for variations due to geometry, Reynolds number, and inclination of the annulus they have included some empirical parameters in their formulation. Johnson et al. [3] consider the rise velocity and the void fraction (volume of gas trapped inside the fluid) of bubbles released in a column filled with a viscoplastic fluid. Finally, Johnson and White [4] consider the gas rise velocity of bubbles in a column in which there is already established a pressure driven viscoplastic fluid flow. In particular they compare their results with a similar air-water system and show that despite an increased viscosity the bubbles rise faster in the air-viscoplastic system. These works are primarily experimental in which they examine the variables of interest to petroleum engineers. In addition, the special case of spherical bubbles has been examined in two papers [5] [6]. In Bhavaraju et al. [5], the authors examine the flow of a single spherical bubble in a Bingham fluid. They consider the case of small Bingham numbers and apply perturbation methods. Using the Bingham number as the perturbation parameter, they obtain a first order correction to the Newtonian solution for flow around a spherical bubble. Their solution, however, is only valid away from the equator of the bubble (great circle in the plane perpendicular to the bubble motion). In addition their solution requires that as the Bingham number goes to zero we recover the Newtonian solution. Whether this is actually the case is unclear. (For example, Beris, et al. [7] claim that this is not the case for a sphere in an unbounded Bingham fluid.) From their 3 solution they calculate the drag coefficient and the shape of the yield surface. In [6] the authors show that spherical bubbles that are initially trapped in a Bingham fluid can be made to rise by applying an oscillating external pressure to the liquid. They use perturbation methods to derive an analytic solution that consists of a flow due to the pulsation and a flow due to the bubble rising under gravity. They have also conducted experiments that illustrate the phenomenon and that support their analytic results. A related problem that has received more attention is the problem of slow flow past a solid sphere in viscoplastic fluids. A key paper is that by Beris, et al. [7]. They examine the fall of a sphere due to gravity in an infinite Bingham fluid. Using a finite element numerical simulation they obtain the flow field and the yield surfaces in the flow. In addition, using asymptotic analysis they obtain analytic values for the drag coefficient in the limiting cases of low and high (near the critical value where the sphere ceases to move) Bingham numbers. Adding to the work of Beris, et al., Atapattu, et al., [8], conducted experiments with spheres in viscoplastics to determine the drag coefficient, the terminal velocity, and the shape of the yielded region. Their results, when compared to those of Beris, et al., are qualitatively the same though some discrepancies did exist. One point of note is that Carbopol, the viscoplastic fluid used by Atapattu, et al., (and the same fluid used in this study) is not a perfect Bingham fluid. It is much more accurately modelled using the Herschel-Bulkley constitutive model. Spheres falling in tubes filled with a Bingham fluid are examined by Blackery and Mitsoulis [9]. They conduct numerical simulations to obtain the flow field for different values of the ratio of sphere radius to column radius. While the shape of the yielded region is qualitatively the same, it does not agree entirely with either that of Beris, et al. or of Atapattu, et al. (In the case of the largest diameter cylinder the walls are sufficiently far away that the result should be identical to the flow in an infinite liquid.) In addition, some work has been done to generalize the equation for the Newtonian drag coefficient for a sphere to one for a sphere in viscoplastic fluids [10] [11]. The results from these papers correspond well with the results of Beris et al. Some work has even been done with 4 mixed slip-stick boundary conditions on the sphere [ 1 2 ] and also to experimentally determine the drag coefficient for other shapes (such as discs, cylinders, cubes, and cones) in viscoplastic fluids [ 1 3 ] . And finally, Li and Renardy have numerically examined deformations of a viscous fluid frop in a Bingham fluid [14]. The experimental work on creeping flow around a sphere in viscoplastic fluids has concentrated more on corroborating theoretical results. Unfortunately insufficient attention has be given to ensuring that the experiments are repeatable and that the results obtained are reliable. Gheissary and van den Brule [15] have shown that in practice viscoplastic fluids have a relaxation time. In general, the terminal velocity of a sphere increases with subsequent releases of spheres. Only if a sufficiently long amount of time between spheres is given (or if the liquid is stirred) does the subsequent sphere fall at the same velocity as the first. In Carbopol, a viscoplastic commonly used in experiments, the result of releasing subsequent spheres is even more surprising. Subsequent spheres dropped in relatively quick succession have an increasing terminal velocity. If a sufficient amount of time elapses (on the order of tens of hours) then a sphere will again fall with a velocity the same as the first sphere, however, if an insufficient amount of time elapses, say only a few hours, the velocity of the next sphere will jump to a value higher than that of all the previous spheres. Thus in any experiment it is extremely important to ensure that before each test the fluid is "reset" to some initial state, thus avoiding the complication of the time dependency. Another closely related problem is that of the propagation of a bubble in non-viscoplastic fluids. Bubbles in Newtonian fluids have been studied for quite some time and there are several books that cover the collection of related work [16] [ 1 7 ] . Also of note, Bretherton [ 1 8 ] has also examined the flow, including effects of surface tension, of long thin bubbles moving in tubes filled with a Newtonian fluid. Moving to the realm of non-Newtonian fluids, Bhavaraju, et al. [5] using their perturbation methods determined the flow around a spherical bubble in a power law fluid. Again their result is only valid away from the equator of the sphere. Vasil'chenko [ 1 9 ] experimentally compares the rise velocities of bubbles in a tube filled with a viscoplastic-elastic fluid against the rise velocities of bubbles in tubes filled with pure viscoplastic and viscoelastic 5 fluids. The former case exhibits an unusual behaviour not seen in the pure fluids. In the viscoplatic-elastic fluid, for large bubbles that fill the tube, the rise velocity is an increasing function of the bubble length, whereas for the pure liquids the rise velocity is constant. Finally, some of the methods used in this study, in particular the variational methods, have already been used to examine viscoplastic flows; however, as yet flows with bubbles have not been considered in this way. The first instance of a variational minimum and maximum principle for Bingham fluids is attributed to Prager [20]. In a series of papers, [21], [22], and [23], Mosolov and Miasnikov used Prager's minimum and maximum principles to study two-dimensional flow in pipes of arbitrary cross section and two-dimensional flow around an arbitrary solid body in a channel. They determined some of the characteristics of the flows and developed ways to compare flows. More recently Glowinski has used the variational methods with numerical calculations [24] [25]. Huilgol has used the variational principles in cases where there is slip present at walls [26], and when inertial effects are included [27] [28]. Frigaard, et al. have used variational methods to study exchange flows and displacements flows of two viscoplastic fluids [29] [30]. We mention these papers here since, in addition to the general methodology, some specific techniques are taken from these. 1.3 Outline of Thesis In this thesis we consider the propagation of a single air bubble in a cylindrical column filled with a viscoplastic fluid. We focus our attention on slow moving and stopped bubbles, and the conditions required for a bubble to become trapped. In chapter 2 we describe the physical problem and write out the full equations that govern the flow. We also nondimensionalize the problem and introduce the important dimensionless parameters. In chapter 3 we simplify the full problem and develop two variational results: a variational inequality that leads to a rate of strain minimization, and a stress maximization result. In chapter 4 we use our variational results to obtain stopping conditions, under which a bubble will not move. The first condition is a general result that applies to axisymmetric bubbles of arbitrary shape. The second condition allows us to use existing solutions of simpler problems to obtain extra stopping conditions. 6 F i n a l l y the t h i r d stopping condition applies to the special case of long c y l i n d r i c a l bubbles. In chapter 5 we derive several results relating the bubble propagation velocity to the physical parameters of the p r o b l e m . F i n a l l y , i n chapter 6 we describe a n experimental setup which we use to observe the the propagation of bubbles i n a c y l i n d r i c a l c o l u m n . W e present our experimental results a n d analyze the data. In addition, we compare the experimental results to our previous analytic results. 7 Chapter 2 Problem Description 2.1 Physical Setup The problem we are considering is that of the flow (or lack thereof) of a single gas bubble in a cylindrical column filled with a viscoplastic fluid. The basic physical setup is shown in figure 2.1. Our entire domain is divided into two connected regions, fi, comprising the liquid viscoplastic fluid bubble interface walls Figure 2.1: Physical setup of the problem. region, and fib, comprising the gas domain (i.e., the bubble). The boundary of fi, consisting of the cylinder walls and the bubble surface, will be referred to as <9fi. The bubble surface on its own will be referred to as <9fi&, and the walls of the cylinder on their own will be referred to as 8 dQ w or simply as the "cylinder walls". Quantities denoted with a hat " ~ " are dimensional while quantities without a hat are dimensionless. We also adopt the Einstien summation convention for the indices. The equations of motion in the liquid region, Q, are „ Dui dp df (u) ttij P= (2.2) 0, OXi is the material derivative, u = (ui, 112, ^3) is the velocity, p is the pressure, where -j^ = pt is the liquid density, is the deviatoric stress tensor for the liquid, and T^J(U) gi = (0,0, —g) is the gravitational acceleration. In the gas region, f^, the equations of motion are „ Dili p where p g 9 is the gas density, and 1 5 dp - M = t _ + df (u) gtij ^ p l - - d J - - ' + ( 2 " 3 ) is the deviatoric stress tensor for the gas. f ij 9i The second invariants of the rate of strain tensor and the deviatoric stress tensor are defined as '1 ^f(u) = \/^ij{u)^ij{u), Tk{u) = \/ ^h,ij(u)f ,ij(u), (2.5) k = £,g k (2.6) where ^ U) = w ^ dxii + - du (2 3 7) is the rate of strain tensor. The liquid, being viscoplastic, is modelled as a Herschel-Bulkley fluid. The constitutive equations for a Herschel-Bulkley fluid are 7i (u)=0 T ,ij{u) t = (2.8) i f T < ( u ) < f y i fA ^ ( « ) V n _ 1 + T777 J 7 W / 7ij(t*) ^ T {u) t > ry (2.9) where fit is the consistency, n is the power law index, and fy is the yield stress. All three of these parameters are positive constants. The gas, considered to be Newtonian, has the constitutive equations Tg,ij{u) (2.10) = fl,gjij(u), where fi is the viscosity of the gas. g Note that for a Herschel-Bulkley fluid the domain fi will be divided into two (not necessarily connected) regions, determined by equations (2.8)-(2.9). One region is where the stress in the fluid does not exceed the yield stress and equation (2.8) is valid. This is referred to as the unyielded region. The other region is where the stress in the fluid exceeds .the yield stress, and equation (2.9) is valid. This region is called the yielded region. The yielded and unyielded regions are separated by a yield surface at which the stress is equal to the yield stress. In addition, we have boundary conditions that must be satisfied on the cylinder walls and on the bubble surface. On the walls of the cylinder we have a no-slip condition: Ui = 0, on dQ . w (2.11) Physically we actually have a free surface at the top of the cylinder; the free surface can rise to account for any expansion of the bubble due to the change in hydrostatic pressure as it rises. However, later on, in section 3.2, we show that in the case of slow flow the bubble can effectively be considered to be incompressible. With incompressibility we no longer need the free surface to be able to rise, i.e., the average height of the liquid in the column will be constant. Moreover, for any finite yield stress and a sufficiently long column, there will be a finite length below the free surface for which the fluid is not yielded. As we will show in Appendix A it is possible to introduce a wall at any point in an unyielded region without affecting the flow. Thus considering a column with a rigid top that is completely filled with a Herschel-Bulkley fluid is a mathematically equivalent problem. On the bubble surface we have a set of jump conditions involving the velocity and the traction. 10 The traction vector is defined as °k,n i{u) = ak,ij{u)n j, b k = £,g b (2.12) where vk,ij(u) = -pSij + T ,ij(u), k = £,g k and Tibj (2.13) is the outward unit normal of the bubble (i.e., pointing into the liquid). The velocity and the tangential components of the traction must be continuous across the bubble surface: ue,i-u ,i = 0, (2.14) °g,n i(u)] h,i = 0, (2.15) [ve,n i(u) - Vg,n i(u)] t ,i = 0, (2.16) g [&e,n i(u) b _ b b b 2 on dfl , where tij and t2,i are any two linearly independent unit tangent vectors of the bubble 0 surface. For the normal component of traction fo,n i(«) 6 - °g,n i(u)] n b 6 | i =| + j-^j , (2.17) on <9fi(,, where n^j is again the outward unit normal of the bubble surface, £ is the surface tension, and (-i- + J - ) is twice the mean curvature of the bubble surface (Ri and 7^2 are the V. Ri R2 J radii of curvature in the principle directions). Finally, having determined the velocity, u, the bubble interface, <9f2(,, whose location is denoted by F(x,t) = 0, evolves according to the kinematic condition dF dF +^ = 0. dt oxi ^ (2.18) Thus for a given bubble, the solution velocity field must satisfy equations (2.1)—(2.18). 2.2 Non-dimensionalized Equations To non-dimensionalize (2.1)—(2.18) we take relevant dimensional scales from the physical problem. For the length scale there are two natural choices, either the column radius or the effective 11 bubble radius. Here we choose the effective bubble radius as the desired length scale. The effective bubble radius plays a more significant role in influencing the flow. This will become evident in the other dimensional scales that we use. Thus we take our length scale to be &i ~ R, (2.19) where R = %J^V and VJ> is the bubble volume at the injection pressure. Another advantage 0 of this scaling for length is that in the experiments we control the volume (and thereby the effective radius) of the injected bubble. The velocity scale is chosen through an approximate force balance. As we are primarily interested in slow flows and cases where bubbles are static or "stopped", there should be an approximate balance between the buoyancy force on the bubble and the viscous force on the bubble. The buoyancy force ~ (pg — p* )gR? and the viscous g force ~ • Since the gas density is not constant, we have written Pa = P* Pg, (2-20) g where p* is the gas density at the injection pressure and p is the non-dimensional gas density. g Since the gas density is very small compared to the liquid density (i.e., pt — p* « pi) we take the velocity scale as * „ t f It turns out that using pi instead of pt~P*g i ( W * ^ \ = sa (2.2!) better choice for U. Finally we scale the pressure hydrostatically as p ~ ptgR. (2.22) We define the dimensionless length, velocity, time, and pressure (xi, Ui, t, and p respectively) by Xi = Rxu (2.23) Ui = Uui = i = ~t, U Ui, V ^ and P = PegRp- / (2-24) (2.25) (2-26) 12 The rate of strain tensor becomes iij{u) = fijr + ^ w h e r = ?7ii( )w e ( 2 2 7 ) The dimensionless constitutive equations for a Herschel-Bulkleyfluidare 7Ji(«) = 0 if T,(U) < B , Ti,ij(u) = (jiu^- 1 + JL^j %(u) (2.28) if r (u) > B, (2.29) and (2.30) t where n,ij{u) = ^—Te (u), R tij n B =^ =^ r (2.31) faU fagR 71 is the (dimensionless) Bingham number. The second invariants of the dimensionless rate of strain and stress tensors are i( ) = \l niij( )iij( ), u u u a n (- ) d 2 re{u) = ]J~T (u)T (u). e>ij 32 (2.33) e!ij The constitutive equations for the gas simply become Tg,ij{u) = Sjij(u) (2.34) where f <W"(«) = ^r 9,ij( )' T 2 n and 8= (- ) U R A g f i im—1 " . 35 (2.36) The non-dimensional versions of (2.1) and (2.2), valid in the liquid region, are F r ^ =-^-6 ^ = 0, OXi 13 -^, d l 3 + (2.37) (2.38) and the non-dimensional versions of (2.3) and (2.4) are e p ° F r -W = -dx-dpg d - f + ^ L e ^ p > dt + ( 2 = 0, - 3 9 ) (2.40) where the dimensionless parameters are Fr = -j=, (2.41) e=$-, Pi • (2.42) figR^ = and 5ij is the Kronecker delta. The boundary conditions (2.11) and (2.14) remain the same: u = 0 on dQ , % (2.44) w ue,i — g,i = 0 across dQ, . u (2.45) 0 The dimensionless traction conditions are n,n ti b n,n t b ~W +Pg+ U,n n b 7" ,n ti 9 6 ~ g,n t T 2 b 2 ~ T ,n n b g b = 0, (2.46) = 0, (2.47) = b P (J^ + ^ , (2.48) across <9f2(,, where B = (- ) 2 pe.gR 49 1 is the dimensionless surface tension, and R\ = and i? = it are the dimensionless radii of 2 It curvature. Note that T jj is in general 0(5) smaller than r^-; see (2.28)-(2.30), (2.34), and g (2.35). Finally, the kinematic condition for the bubble surface, (2.18), remains the same 14 Chapter 3 Variational Formulation 3.1 Introduction The existence of the yield stress, and hence the possibility of a yield surface, makes the classical approach of directly solving the Naiver-Stokes equations impractical. The location of a yield surface is determined by the stress field and the yield surface in turn determines the boundary conditions for the system. Furthermore, the yield surface need not be static, nor is it a material surface. In essence we have a complex free boundary problem, for which an analytic solution can be obtained only in special cases. Another approach which can be used is to formulate the problem as a variational problem, where the actual flow field distinguishes itself from all other possible flows in that the actual flow minimizes or maximizes some functional. While the variational formulation does not result in an explicit solution for the flow field, we can nevertheless learn a lot about the behaviour of the solution. The main motivation for using a variational formulation is that it eliminates the need to know the location of the yield surface beforehand. 3.2 Simplification of the Full Problem For convenience we rewrite the full non-dimensional problem. In the liquid region, the nondimensionalized equations of motion are Fr Dt — dxi — 6 dxi (3.1) (3.2) 15 and in the gas region the equations of motion are 2 i -m Du e p ° F r P -d - d = x ^ ^ p to o\ - , jr^TijC") x e i -dx—> + 5 ( 3 ^(P^)=0, + where F r = — 7 = is the Froude number, B = = - i i ^ - is the Bingham number. 3 ) (3.4) is the Kronecker delta, e = -f, 8 — ^-r — , and T r For a Herschel-Bulkley fluid the (dimensionless) constitutive equations are jij(u) = 0 ri,ij(u) = (iiu)"- + 1 %(u) XTI{U)<B, (3.5) if ( u ) > 5 , (3.6) T / and for the gas the (dimensionless) constitutive equations are T ,ij{u) = 6iij(u), (3.7) g where '1 7 ( « ) = \l ipij ( )%( ), u u a n d Te{u) = y^T ij(u)Te,ij{u). (-) 3 8 (3.9) ei As for boundary conditions, we still have zero flow at the cylinder walls: Ui = 0 on dQ . (3.10) w Across the bubble interface, <9f4, the dimensionless velocity and traction conditions are ue,i - u g>i = 0, (3.11) Te,n t! - Tg,n ti = 0, b (3-12) b e,n t - g,n t = 0, T (3.13) T b 2 b 2 -Pi +Pg+ n,n n ~ Tg,n n = P b b b 16 b + j^J , (3.14) where n& is the outward normal to the bubble surface, t\ and ti are two linearly independent tangent vectors of the bubble surface, and R\ and R2 are the dimensionless radii of curvature in the principle directions. And finally, the kinematic condition for the bubble surface is dF ^ Ot OF + u " i l = 0, (3.15) OXi where F(x, t) = 0 is the location of the bubble surface. For slow flows we can ignore the inertial terms. The inertial terms scale as ~ u f , while the remaining terms scale as ~ U{. Since we are interested in the phenomena of stopped bubbles and "nearly" stopped bubbles, the inertial terms can be considered to be an order of magnitude smaller that the other terms. Also for our system of polymer solution and air S ~ 1 0 - 6 and e ~ I O . Thus, as a first approximation, we ignore the inertial terms and the terms containing -3 S and e. Equations (3.1) and (3.3) become liquid: 0=-- d H gas: ^ i3 OXi OXj , (3.16) 0= - ^ - (3.17) OXi respectively. Equation (3.17) implies that the pressure in the bubble is constant (really p = p(t)). The dimensional version of equation (3.4) can be written as dt id p g P g dp oxi y di_\du V dt dxi J + 1 dp dp g p g 1 2 g g 1 = o dx. | dp dt pc where c — y dxi fdp duj _ Q dxi dp diii di dii = 0, (3.18) is the speed of sound in the gas. The main change in the pressure is due to the static pressure of the liquid, as the bubble rises; this changes like ptgU• Examining the relative sizes of the two terms in (3.18) we have I f ~ 1 1 dp PgC di U R dibj dxi ect 2 g gR_ and ^"WH ~ (10)(10- ) _ 2 17 3 So Thus we can make the approximation of incompressibility. dui „ =0 OXi dui or p = 0, OXi (3.19) in the bubble. Lastly the traction conditions at the bubble interface, ri,n = btl dflt, simplify to 0, n,n t = 0) (3.20) a n ( - ) d 3 b 2 -Pt +Pg+ since T jj g is 0(6) n,n n b = b + -^j S 21 (3.22) smaller than r ^ . Equations (3.16), (3.17), (3.2), and (3.19), with the constitutive equations (3.5)-(3.7), and boundary conditions (3.10), (3.11), and (3.20)-(3.22) comprise the simplified system for which we will develop a variational formulation. 3.3 Variational Inequality (Rate of Strain Minimization) For the variational inequality we obtain a result where the solution function (in this case a velocity field), chosen from a given function space, satisfies certain criteria with respect to the other functions in the function space. The space of functions we consider, denoted V, is the space of all vector-valued functions v = (vi,V2,Vs) Vi ~ oxi such that G C°°(fi), =0 (3.23) infi,and Vi = 0 on dn . w (3.24) (3.25) It is clear that the actual solution u & V. 3.3.1 A Preliminary Inequality First we will derive an inequality which will be employed in our variational derivation. Let itj be a velocity field that satisfies equations (3.2), (3.5)-(3.11), (3.16), (3.17), and (3.19)-(3.22), 18 and let vi be any admissible velocity field from our functional space, V. In the unyielded regions, equation (3.5) holds and the quantity -T {u)%{u) = 0 = ^ 7 ( « ) l n _ 1 itij ^T£ JJ(U)7JJ(W) ) satisfies 7 y ( « ) 7 y ( « ) + B^(u), since both jij{u) = 0 and 7(14) = 0. Also the quantity ^T£ ij(u)jij(v) (3.26) satisfies : ^ ^ , y ( « ) 7 i i ( « ) < T<(«)7(«) < Bj{v) = ^7(u) n _ 1 7 i i ( « ) 7 i i («) + (3-27) Here we have used the Cauchy-Schwarz inequality otijfiij < ^Jcx-ijdij A/PijPij, and that while we do not know the exact value of Tg(u) in the unyielded region we know that it is less than the yield stress, B. Again jij(u) is simply zero. In the yielded regions, equation (3.6) holds, and the quantity ^T£ ij(u)jij(u) satisfies t \n,ij(u)^ {u) = i (i(u) ~ n x i:j 1 = 27(«) = ^7(«) iij(u)iij(u) 7ij(«)7y(«) + n _ 1 n _ 1 + ^TJ) 7»j(«)7y(«) IB ^^fju)^^^^ (3.28) +Bi(u). And the quantity ^T( ij(u)jij(v) satisfies t \T^ {U)^ {V) 3 13 = ^ (i{u) n 1 = 27(«) n _ 1 + ^-^j 7 i j ( « ) 7 i i ( « ) 1 1 _z? 7ij( )7ii(») + 2^)^ '^ ^* ^ ^ u ij < liiuT-^iu^jW u + Biiv), J v (3.29) where we have again use the Cauchy-Schwarz inequality in the last line. Comparing equations (3.26) and (3.28), we can see that the equality ^^,y(u) „(u) = )p(u) - iij{u)ii (u) n l 3 7 19 + Bj{u) (3.30) holds in both the unyielded and yielded regions. Similarly, from the inequalities (3.27) and (3.29), we have that the inequality (3.31) also holds in both the unyielded and yielded regions. Subtracting (3.30) from (3.31) we have that 2W u ) (7y(«) - 7 « ( « ) ) < \n,ij(u)iij(v 7J7H" - u) < ^ 7 ( w ) n _ 1 - 7ij(«)) + Bi(v) - Bi{u) 7 v ( « ) 7 « ( « ~u) + Bj(v) - B-y(u), (3.32) since jij(') is linear in its argument. Inequality (3.32) is crucial to the variational formulation. Inequality (3.32) avoids the problem of determining where the yielded and unyielded regions are located - it is valid throughout the liquid domain fi. Later in our variational formulation we will return to this inequality. 3.3.2 Derivation of the Variational Inequality Again let m be a velocity field that satisfies equations (3.2), (3.5)—(3.11), (3.16), (3.17), and (3.19)-(3.22), and let Vi be any admissible velocity field from our functional space, V. Multiplying equation (3.16) by (UJ — ui) and summing over the index i we have ^— + S 1 (Vi - Ui) + — [Te,ij(u)(Vi - Ui)] - T( j(u) — (Vi - Ui) i3 ti where we have used the chain rule to rewrite the last term. 20 (3.33) Since T^J(tt) is a symmetric tensor the last term of equation (3.33) can be written d e,ij(u)-^-{vi T 1 - d = -T^(«) —(^ dxi 1 - d + -^n,ij{u) — (vi - u^ dxi = e,ij(u) T 2W ^- ^ U dx-^- + 3 U ^) 1 = 2^j()7ii(« - «)• T Furthermore, due to incompressibility (3.34) u d(vj-Uj)u d x = 0, so the first two terms of equation (3.33) can be written - + <5i3^ (Vi - Ui) + [Ti j(u)(Vi ti - Ui)] dp d dx (Vi - Ui) - {v - Uz) + — d Q \p{Vi - Ui)} - (v - U ) + [T (u){Vi - Uj)]. d [n {u){Vi - Ui)] 3 - x 3 e>ij 3 tij d = -7^- d \pSij(Vi - Ui)] - (v d U) + \[-p5ij 3 = -(V 3 — 3 = -(V - U ) + — 3 where a(. ij{u) t dxi [T£,ij(u)(Vi - U) + — 3 + T j(u)] tti [oe (u)(vi 3 (Vi - Ui) (3.35) - Ui)] , tij - Ui)] = —p5ij + T^JJ('U) is the dimensionless traction. Thus equation (3.33) can be be written 1 2 e,ij( )iij( T u v -u) = -(v 3 d u ) + — [o ,ij(u)(vi - u^]. 3 (3.36) e Now we use the inequality (3.32) on the left hand side: 1 d •^i{u) ~ iij{u) f {v n lA A ij -u) + Bj(v) - Bi(u) >-{v 3 u ) + — [a (u)(vi 3 etij - u^]. (3.37) Finally we integrate equation (3.37) over the domain Q to obtain a(u,v - u) + j(v) - j(u) > - [ (v -u )dQ+ Jn 3 3 21 ( -^—[(j ij{u){vi-Ui)]dVL Jn OXJ lt (3.38) where a(u,v) = \ [ j{u) -yi (u)ji (v)dfl, n and 1 j j j(u) = B f ( « ) 7 da Applying the divergence theorem to the last term of the right hand side of equation (3.38) we have a(u,v - u) + j(v) - j(u) >- (v -u )dQ 3 + / 3 ae,ij{u)(vi - Uijrij ds, (3.39) JdQ. Jci where rij is the outward unit normal of dft. Splitting the boundary into the wall boundary and the bubble surface, and noting that rij = —ribj on dQb, we can write a(u, v -u)+ j(v) - j(u) > - {v - u3) dtl 3 Ju + ae j(u)(vi - }i Ui)rij JdU ds - w a (u)(vi - eiij Ui)n j b JdCl ds (3.40) b Recall that n j is the outward unit normal of the bubble. Now on dft b the veocities, Uj and w Ui are defined to be zero. Thus the integral over 8Q in equation (3.40) vanishes. Futhermore, W on the bubble surface, dUb, the traction must satisfy conditions (3.20)-(3.22). We can write / &e,ij(u)(vi - Ui)n ds Jdn = b>j b / [a^ n {u){v Jdn nb b nb - u) nb b + o W 1 ( t t ) K 1 - u ) + a£ n (u)(v h t bt2 - u )] ds. t2 t2 (3.41) From (3.20) and (3.21) the last two terms are zero, and for the remaining term, using (3.22) we can write / &e,ij(u)(vi - Ui)n j JdU b b ds= f Jdn -p L g b + (3 (-J- + J rC J 2 22 [v nb - u ) ds. nb (3.42) Now, since the velocity is zero along the cylinder walls, / -Pg(v nb ~ U n ) ds = / p (vi b Jd<n g b / - / Jn - Ui)ni ds w p {vi g - Ui)(rii) ds " \Pg( i - i)] V = Pg g Jan b = p (vi - Ui) {-n ,i) ds + / Jon b v Jn —0 u V • (t»i - d£l Ui) (from incompressibility; see (3.2)) Thus the integral over the bubble surface, equation (3.42), reduces to / atij{u)(vi Jan ds = [ @ ( + Jan \ i - Ui)n bjj j (VJ Ri) K b b Ui)n j b ds (3.43) Finally, from equations (3.39) and (3.43), we arrive at the final variational inequality for the problem: a(u,v - u) + j(v) - j(u) > L(v - u) Vu G V (3.44) where a(u,v) = \f ^{v.) - %{u)^j{v)dn, Jn n l 1 j(u) = B / 7(u)dft, Jn L[u) = - I u dttJn 3 and [ p(i- +^-\ Jan \ i 2j K K (3.45) (3.46) U i n b i ds. (3.47) b It should be noted that a(-, •) is linear in its second argument and L(-) is linear, while j(-) is non-linear. In some of the literature a(-, •) is referred to as the viscous dissipation rate and j(-) is called the yield stress dissipation rate. 3.3.3 A n Alternate Formulation Equation (3.44) requires us to find a u G V such that for all other v G V the inequality (3.44) holds. Alternately we can formulate the problem as a pure minimization problem. Consider the function H(u) = -L- (u,u) a = —±—r 23 / T H " - 1 ^ ^ ' ^ ) ^ . (3-48) The Gateaux derivative of H in the direction v is 6H(u;v) = \f T H ^ H J H ^ H ^ Jn 1 = a(u,v). (3.49) (3.50) Furthermore H is convex. Since H is Gateaux differentiable and convex we have from [25] that H(v) - H{u) > 8H(u; v - u) -a(v,v) n +1 ' v -a(u,u) > a(u,v — u) n+1 (3.51) Substituting (3.51) into (3.44) and recalling that L(-) is linear we obtain 1 n + -a(v,v) 1 v \-ra(u, u) + j(v) - j{u) > L{v) - L(u) n+ 1 —^—a(v, v) + j(v) - L(v) > —!—a{u, u) + j(u) - L{u) n+1 n+1 Vu G V. (3.52) So the true velocity field, u, distinguishes itself from all other velocity fields, v € V, in that u minimizes the functional J(v) = -^—a(v, v) + j(v) - L(v) n + 1 (3.53) over the functional space V. Also in equation (3.44), if we let v = 2u then a(u,u) + j{u) > L(u), (3.54) where from equations (3.46) and (3.8) we have written j(2u) — 2j(u). Taking v = 0 in equation (3.44) we obtain —a(u,u) — j(u) > —L(u) =>• a(u,u) + j(u) < L(u). (3.55) Equation (3.54) and (3.55) imply that the true velocity field, u, must satisfy a(u, u) + j(u) = L{u). (3.56) Thus, from (3.56), we also know the value of the functional J at the minimum: J(u) = TI a(u,u). n + 1 24 (3.57) 3.3.4 Existence/Uniqueness For the special case ofra= 1 in the Herschel-Bulkley constitutive equations case of a Bingham fluid), the variational equality (3.5)-(3.6) and the minimization problem (3.44) (the (3.53) have a unique solution in the functional space. This is shown by directly applying Theorem 4 . 1 and Lemma 4 . 1 from chapter 1 of Glowinski [ 2 5 ] . For the case of n ^ 1, to the best of our knowledge, no general existence/uniqueness results exist. However, if existence can be shown, it is then easy to show uniqueness using the method used by Glowinski. Throughout we will assume that there exists a unique solution to 3.4 and (3.44) (3.53). Stress Maximization Principle The earlier result, (3.53), derived from the variational inequality is also known as a rate of strain minimization. Similarly there is another formulation that leads to a stress maximization principle. Consider a stress field ot^j. We say that the stress field = —p~e8ij + or equivalently the pair (p~g, f^,-) is admissible if n _ U — -T; dpi OXi h,n ti = 0 h,n t = 0 b b -PI + h,n n 2 b b = -P9 dfi j ti r + -T; OXj Oj f (l^ 3 + 1 in Q, (3.58) 3 - \ on dQ , b (3.59) on dn , (3.60) on dn . (3.61) b b R Let T denote the set of all admissible (pe^e^j). T h e o r e m : If U{ is the actual solution of the classical problem, then the true stress field, r ^ - , maximizes the functional F{h,ij) = r ^ r - ^ - r / (\re -B\ + f t 2 n + i n + 1 B)l +l Jn dn+ f f^7ij(u) dQ, V {p , f ) t tjij G T. Jn (3.62) 25 Proof: T h e p r o o f is taken from [30], s u i t a b l y modified. F(T ) - F{T ) ltij = I ttij } [(\t -B\+ti-B)^ n + i -(\T -B\ + 1 l + T i - / t - B ) l - n,ij)iij(u)d& Jn + l (3.63) C o n s i d e r the integrand 1 2^+1 n n + (\f t 1 -B\ + f - B)n - (\ +1 t + T - B ) - -B\ rt + - {f 1 t ltij - Ttjj)iij{u). W e w i l l show that at any point i n fl, I > 0, a n d thus the integral (3.63) is greater t h a n or equal to zero. T h e constitutive equations for a H e r s c h e l - B u l k l e y fluid c a n also be w r i t t e n as iii = 0 iir <B (3.64) ii (3.65) e > B TI Case 1: 'jij = 0 =>- T£ < B Then 1 I= Case 2: 7 > 0 =4> r n \f e -B\+f - B ) n t + l > 0. > B t T h e n using (3.65), 1= n 2^+1 n + 1 (|f 1 -B\ + B ) l TI- + l - (\T -B\ E + t -( Case 2a: B ) ^ r - T e + 1 -B)^(f i t i j -n } i j )^-. (3.66) fg>B Since we also have > B, (3.66) simplifies to n n + 1 [f - B)^ +l t - [rt - B)^ ] +1 26 - (r e - B)±(f ttij - T ) tiij re ' (3.67) For all Tijj of a fixed magnitude I is minimized when fg^j \\ Tg^j, i.e., fgjj = Xrg^j. We can also write Tg = 9B ; (3.68) for some 9 > 1. Similarly, = AT i = (3.69) X9B, where A0 > 1. Substituting (3.68), (3.69), and fgjj I > I*(A) = B^ = XT^J into (3.67) we obtain ((A0 - l ) n - (9 - 1 ) ^ ) - 9(9 - 1)*(A - 1) +l +1 +1 Looking for the minimum of I*(A), 9(X9 - riA l)n - 9(9 - 1)» when A = 1. Also |£-I f l i V(tf-l)i->0 + since A0 > 1. Thus /*(A) is minimized at A = 1. Therefore I > I*(l) =0. Case 2b: fp<B Also TI > B, so (3.66) simplifies to I = -- n - ( 7 7 - - (T* - B)n(f eiij - Again for all fg^j of a fixed magnitude I is minimized when fg^j Ti ) tij T£,ij (3.70) || T ^ J J , i.e., f ^ - = A T ^ J . We can again write TI = 9B (3.71) for some 9 > 1. Similarly, = AT/ = A05, 27 (3.72) where A0 < 1. Substituting (3.71), (3.72), and f^ij = Xrg^j into (3.70) we obtain I>f(A)-B» + 1 - 0 ( 0 - 1)~(A- 1) — — ( 0 - l)n [n + 1 +l 7* (A) is linear in A and decreases as A increases. The maximum allowable value of A is restricted by A0 < 1 => A < \ 0 Thus /*(A) is minimized at A = ^. Therefore >0, since 0 > 1. • The above proof does not make use of any the conditions of an admissible stress field. Using these conditions, we can obtain a useful relation. The derivation of the relation is almost exactly the same as the derivation at the beginning of section 3.3.2. The admissible stress field satisfies (3.58) and the fluid is incompressible (i.e., ^ dpi 0= = 0), so for the true velocity field Ui, we have df ij it Ui - U 3 + Ui — dxi dxj Since f^j is symmetric (see section 3.3.2, equation (3.34) ), Combining (3.73) and (3.74), we have 1 d 2^,u7»j(w) = —u + g^r;[(-Pt ij + t,ij) i\s f u 3 Integrating this equation over the domain fi, and applying the divergence theorem to the last term, we obtain T; f n,ijiij( ) Jn u 1 = - [ u cinJn 3 28 [ {-peSij+ fe )uin jds. Jdn yij b b (3.75) The sign in front of the last term changes since we have changed from the outward normal of Q on the bubble interface to the outward normal of the bubble; n ^ = — n* on dQ , where rii b b is the outward normal of Q. On the bubble surface, we consider a local coordinate system, {rib,ti,t2}, with one axis, (rib), perpendicular to the surface of the bubble. Then / {-Pthj + n,ij)uin j JdQ ds = b (-peS + inb JdQ b h )uids jinb b = / [(-Pt + h,n n )u b b + f u nb ittinb + f tl u] £Mnb ds. t2 (3.76) Prom the conditions for an admissible stress field (3.59), (3.60), and (3.61) we can write / (~PAJ + f j)uin j Jon iti ds = bt b (3.77) Uiiih^ds Jou b Now since the velocity is zero along the cylinder walls and the pressure in the bubble is constant (see (3.17)), / -p Uin ds= JdQ g / PgUi(-n^i) ds + / PgUiUids JdCi JdCl bji b b = w / PgUiUids JdU dU (p Ui)dQ g = 0. (from incompressibility; see (3.2)). Thus combining this result with (3.77), (3.75) becomes 7i z f Jn h,ijii3(u) dQ, = - / u dQ- Jn [ 3 B Jan h ( + \i K 2j K (3.78) Uin ds. b>i Therefore we can we can write the functional (3.62) as F(f ) eij = * \ [ (\r - B\+f e e B)n + i dQ - 2 f u dQ - 2 / Jo. 29 z ^ (i _+ i _] n \i 2/ U i Jdn K b K b t i ds. (3.79) The last two terms on the right-hand side of (3.79) are constant for a given bubble, only the first term depends on f ^ . Thus the stress minimization principle can be stated as: For a bubble of given volume and fixed shape, the true stress field maximizes the functional G(f ) = - - L - - ! ! - f (\f -B\ + f - B)^ +1 £tlJ e 2 n + i n "T 1 t JQ 30 dfl, V (pi, f ) ttij € T. (3.80) Chapter 4 Stopping Condition Results We are primarily concerned with slow moving and stopped bubbles, and conditions which determine whether a given bubble will propagate in the viscoplatic fluid or will become trapped. In this chapter we use the variational results of the previous chapter to derive conditions under which a bubble will not move. 4.1 4.1.1 First Stopping Condition Condition O n B For N o Motion From our variational formulation we have that the actual velocity field, u, satisfies the equality (3.56). The term a(u,u) = \ J 7 ( w r 7 * j H 7 i j ( ) ^ > 0 _1 u so we have the inequality j(u) < L(u) Thus, in a situation where the actual solution velocity is non-zero B f 7 ( u ) d f i < - [ u dflJn Jn 3 B l [ 8 ( + Jdn V-^i b hn ^{-k In^dQ k b ) ^-2/ + Uin jds, b K l ) ^ d s Or equivalently, for B JQ vs dfi > ~ au^ev I I n P ( i r + -k) * ViUb d f i(v) dfl b J 7(1;) dii Q n 31 ds \ (4.2) the only possible solution is one with u = 0 identically over Cl. These results hold in general for any arbitrary bubble shape and bubble velocity. Using (4.1) we could, for example, determine, for a given Bingham number (and a given volume), shapes that will not move, if such shapes exist. Prom (4.2) we could find, for a given bubble (shape and volume), an estimate of the "critical" Bingham number, above which the bubble will not move. From this point on we will restrict our attention to steadily moving axisymmetric bubbles. 4.1.2 Surface Integral Term Consider an axisymmetric bubble moving along the axis of symmetry in a cylindrical column, as depicted in figure 4.1. The radius of the bubble at a height z is given by the function r = f(z), z =L z r = z = z = Z- z =0 r Figure 4.1: Axisymmetric bubble centred in the column. with f(z) =0 for z > z + and for z < Z-. The second term in (4.1), fni(u)dfl 32 (4.3) involves an integral of the mean curvature of the bubble over the surface of the bubble and is, therefore, explicitly dependent on the shape of the bubble. For the analysis of this term we make use of some basic differential geometry results, taken from [31], to obtain an expression for the mean curvature in terms of the functional shape of the bubble. The calculation of these results are contained in Appendix B. In this derivation we assume that we have an axisymmetric steady bubble travelling along the axis of symmetry of the column. In this case the curvature of the bubble is only a function of z. Let Then, using incompressibility, we can write d . , . . da(z) da(z) (a(z)ui) = Ui = udx dxi dz — 0 (4.4) 3 Integrating (4.4) over a small horizontal slice of fluid, as depicted in figure 4.2, and applying V dV w R c Figure 4.2: Annular slice of volume, V, over which we are integrating. the divergence theorem to the left, hand side we obtain I a(z)uin ds = / u —-—dv, Jdv Jv uz. Vji 3 (4.5) where nv,i is the outward unit normal of the fluid element V. It should be noted that for the portion of dV that lies on the bubble surface, ny,i has opposite direction to n^j, the outward 33 unit normal of the bubble. Expanding the left hand side of (4.5) into the separate sections of of dV we have / ' a(z)uinv,ids + / JdV ' b / a(z)uinv,ids+ JdV+ ' a(z)uiny,ids JdV- f j / a(z)Uinv,ids JdV + f / u Jv = w rz+Az I 2Txa(z)f(z)-Jl Jz r / JdV+ + ( / ' ( z ) ) Uin ,idz+ 2 v da a z z a(z)u ds 3 f , x , / a(z)u ds Jdv- — ( )j —-—dv 3 3 n + 0= f / Jv da{z) , u —-—dv. dz 3 Here we have used the result (B.13) from Appendix B. Dividing by Az and taking the limit as Az -> 0+, 2ira{z)f(z)y/l + {f'(z)) 2 mn i + ^ Vl => 27ra(z)f(z)^l a(z)u ds^J 3 + (f'(z)) u n ,i 2 t v = u ds = -<*{*)§^ (f 3 ) ' U3<is (46) where S is the horizontal surface at height z that passes though the fluid region (seefigure4.3). - r = R c cylinder wall Figure 4.3: Top view of the horizontal cross section S and S . b Now for the entire cross section of the column, S Li Sb, shown in figure 4.3, the flux through this surface must be zero. In fact, across any complete cross section of the column the flux must be 34 zero. Thus / u ds + / u 3 ds = 0 Js Js 3 b} b *i{fs ) U3ds - l { L = U b > d ° ) ' ( 4 ' 7 ) where u ^ is the vertical component of the velocity of the bubble surface at a given point. For b a steady bubble (i.e., for a constant interface shape and constant interface velocity), u 3 is just bt the constant bubble rise velocity which we denote U . Then b / u Js h]3 ds = U b (l)ds Js b b = U x (cross sectional area of bubble). b For an axisymmetric bubble with profile r = f{z), as in figure 4.1, the cross sectional area is TT(/(Z)) . 2 So from (4.7), = -2TrU f(z)f'(z) b Combining (4.8) with (4.6), 2ixU a{z)f{z)f{z) 2wa{z)f(z)y/l + {f'{z))*Uin i = b Vi ^(zKn^4 ' b 1 t a ( z ) / ' ( z ) y/i + (f(z)r (4.9) Using the result (B.12) from Appendix B, for the mean curvature of an axisymmetric surface, t \ a(z) R( = Q M M 1 [-5- + \R -5- off-if') = -P— RJ X -I 2 - 3 - ( A s m ( 4 - 1 0 ) / ( ( / ' ) + l) = 2 2 where / = f(z). Also since there is no net flow out of the cylinder / u dQ = — u 3 dQ Jn Jn 3 bt b = -U [ (1) dQ Jn b b = -u v b b =>U = -^r [ U dQ b 3 n Jn 35 (4.11) where V is the volume of the bubble. b From (4.9), (4.10), and (4.11), f p(^r J \Ri + ^r)u n ds RJ "' l P anb bi = -^-( [ u dSl) I V \J JJ '^" f 33 l 2 b n ~^ f dnb 7 + l) /((/')2 ds X ) 2 (4.12) From Appendix B the surface integral over the bubble is given by (B.13), r2n / pz+ ds= ds= /r JdQh ionb r JO fyjwy+idzdd. (4.i3) Jz- Therefore (4.12) simplifies to I e(±. ±)u ^ + Jan b 4.1.3 \Ri in „,*,) r/w-(/•)',-'>* = s V R2J \J b J J_ n ((y/p z + (414) )2 1 Bubbles Which Never Move Substituting the expression (4.14) into (4.1) we obtain the result that for a non-zero solution to exist D B j u dn /n7W^ n 3 / J u dQ \ r*+ « 27T/3 | N 3 f V v \f -y(u)dn) J _ b n f _ ((y/)2 z + ( / 1 0 2 _ 1 } ^ ) § <-Jt^» f..^ r n n - ( / T - D J j {u)dn v j_ y al b (^/)2 z + 1 )2 ( 4 , 5 ) y For a stationary or upward moving bubble, the net flux of fluid through any horizontal surface must be less than or equal to zero. So if the bubble is actually moving upwards the integral of the vertical velocity will be strictly negative (and equal to —U V ). Also 7(11) is always greater b b than or equal to zero. And for a bubble that is moving (i.e., that has yielded the fluid in some region) the integral of 7(u) 7(11) over ove fi will be strictly positive. Therefore the term — m (4.15) is strictly positive. So if 1 2^ ' p .nr/ - { f f - » «/') + l) 2 d z < 0 ( 4 1 6 ) 5 then it is impossible for (4.15), and hence (4.1) to be true, and thus the only possible solution is one which is identically zero. 36 The expression in (4.16) can, in theory, be made to be less than or equal to zero if the shape of the bubble is such that the integral s = r n n - ( r ? : i ) «/') + i) 2 i z . (4 17) 5 is sufficiently large. It should be noted that mathematically it is always possible to construct a shape for which this term is positive. For a given shape, if S is not positive, then reflecting the shape in a horizontal plane would change the sign of the integral, resulting in a positive value. (Essentially we are applying the transformation z —> — z / —>• / , / ' —>• —/', and/" —>• / " , which we can see changes the sign of S.) To show that in practice it is also possible to' have both negative and positive values of the surface integral term, S, in figure 4.4 we show two actual bubble profiles (bubble velocity is upwards) for which S has opposite sign. f(z) (cm) Figure 4.4: Two bubble profiles. For the solid line S = -0.010, and for the dashed line S = 0.012 (5 is nondimensional). 4 14 Bounding J n We can use the result of the previous section to determine which shapes of bubbles will not propagate for a given yield stress. Another problem of interest is the opposite problem: for 37 what yield stress will a given bubble no longer propagate. For this problem we use the stopping condition (4.2). Including the result (4.14), condition (4.2) becomes i » allv > |^ /-7(tl)dfi ^ ( I- ¥ Vr J _w i 1 n b z n ',-' H> ((J/)2 <««> ) + 1 ) 2 In order to use (4.18) we need to find an upper bound for the term -/n"3<*n ( 4 1 9 ) J j(u)dfl a While we could use our previous result that — f u dfl = U V , we instead relate n f 7 ( t t ) dfl n 3 b b J N u dfl to 3 so that we can obtain a uniform upper bound that does not change as we move to the limit of a stopped bubble. To do this we consider the fluid domain fl divided into three regions: the fluid above z , the fluid below z- and the fluid between Z- and z . + + For the fluid above z + U3 = J ^r(x,y,z)dz. This follows from the fundamental theorem of calculus and the fact that u = 0 (in particular u = 0) at the top wall of the cylinder (i.e., at z = L). So for a horizontal cross section at 3 height z > z + j u dxdy 3 = J (^j ^-(x,y,z)dz j S = ["([ = J ^(x,y,z)dxdy\ ~^(^f dz u (x,y,z)dxdy j S] 3 dxdy (4.20) dz. The integral over x, y is simply the net flux through a horizontal cross section. Since the fluid is incompressible, we must have conservation of volume in the region above this plane. Thus the net flux across this plane must be zero: j u dxdy 3 = 0. Jx,y Similarly for the fluid below z. U3 = -^r{x,y,z)dz. 38 (4.21) For a cross section at a height z < z_, using conservation of mass of the fluid below the surface, j u dxdy = J (^j 3 ^-(x,y,z)dtj dxdy = I\IJii ^ ) ~ {x I Jx,y d U3(x,y,z)dxdy j Ql=(^f = J )dxdy S z dz u dxdy = 0. (4.22) 3 Thus the only contribution to the numerator (4.19) is from the fluid between z and z_. So + u (x,y, z) dx dy\ [ [ [ Jz- \Jf{ I f{z) <x +y <R f u dtt= Jn + 3 3 2 2 2 (4.23) dz. 2 To relate this to the rate of strain we again use the fundamental theorem of calculus to write the velocity in terms of its derivative, and then relate the derivative to the (second invariant of the) rate of strain. Let zi = ~^ ~ be the midpoint of the bubble length and let fi be the z+ z + region of fi above zi andfi_be the region of fi below zi, then for z > zi we have [ Jf{z) <x +y <R 2 2 2 u dxdy = [ I [ ^-(x,y,z)dz\ dx dy Jx,y W z i J 3 2 a z du 3 <-[ Jzi dz (I z) <x +y <R \Jf{z 2 2 2 2 1 lf(z) <x +y <R 2 2 [If 2 (x,y,z) dx dy •y(u(x,y, z)) dx dy 2 2 2 Similarly for z < z\ we have u dx dy < 3 f(z) <x +y <R 2 2 2 2 39 dz V% 2 <— I 1/ -y(u{x,y,z))dxdy / i vv 2^ Jzi \Jf(z) <x +y <R c 1 2 dz dQ. ) dz Hence from (4.23) we obtain ^(z -Z-)J Mu)d(l. + (4.24) a Using (4.24) we get a bound for (4.19), f i{u)dfl ~ 2V2 a Thus from (4.18) we can write Note that this no longer depends on the velocity field. So if (4.26) is true, then the only possible solution is the zero solution. 4.2 Second Stopping Condition In the previous section we used the rate of strain minimization result to obtain the stopping condition. For our second stopping condition we use result (3.80) from section 3.4: The actual stress field, Tg^j, maximizes the functional fj\*l = - over all admissible stress fields (pi,?^) G T. ~ B \ + Tl- B) n l +l dii, Clearly the integrand in (4.27) (4.27) is always non- negative. Thus VfeyeT. G{f )<0 tiij The integrand in (4.27) will have a non-zero contribution to (4.28) G(fg^j) if and only if T > 73 at some point x G fl. Thus if we can find an admissible stress field with G{T ,ij) t 40 = 0, (4.29) then for the actual solution, T^, 0 = G{f ) < G{r ) t>ij ltij < 0, (4.30) and therefore G(TIJJ) = 0 and r < B everywhere in fl. Thus, if for a given problem we can find an admissible stress field then we can obtain a stopping criterion. Possible admissible stress fields can be obtained from the solutions of simpler problems. The solution of a problem with the same geometry but with a Newtonian fluid (or a Power Law fluid) would provide an admissible stress field. We could also extend the fluid domain to be infinite. The stress field solution to this new problem, suitably truncated, will also be an admissible stress field. Moreover we could combine these two simplification and consider a problem with the same bubble geometry but in an infinite Newtonian fluid. As an example we consider the case of a spherical bubble and obtain a new stopping condition based on this method. 4.2.1 S e c o n d S t o p p i n g C o n d i t i o n for a S p h e r i c a l Bubble For the case of a spherical bubble, the admissible stress field that we will consider is the stress field for a Stokes flow solution for a spherical bubble rising under the influence of gravity in an infinite Newtonian fluid. We truncate the stress field appropriately to fit inside our finite domain fl. We denote this truncated stress field for the Newtonian case as TE-. Since the Newtonian stress field is obtained directly from solving the Stokes equations in an infinite domain, rf- will satisfy the momentum equations in fl. Also the conditions (3.59)— (3.61) are imposed boundary conditions in Newtonian problem. Hence rf- is an admissible stress field. Stokes F l o w Solution for a Spherical B u b b l e i n a N e w t o n i a n F l u i d The derivation of this result was first done by Rybczynski, and independently by Hadamard in 1911, and can be found in several standard textbooks [32] [33]. Using a spherical coordinate system (r, 9, <f>) with the variables defined as in figure 4.5, the 41 9 f (r,e,4>) / T / / / / / / Figure 4.5: Spherical coordinates. solution of the non-dimensional problem (using the scalings from section 2.2 with n = 1) is 3r (4.31) 6r For the Newtonian fluid (4.32) so cos 9 2cos6> T N £,rr ,T A f = - - f ut e,r<f> — ) e.r9 — > N N 3r 2 ' T U 3r T N C0S#| V3r cos 9 2 -n (4.33) 2 With our scalings the surface of the bubble is given by r = 1, so T N < -^=, Vxefi. (4.34) ~ V3 Thus with B = ^ G(T£ ij) t we have r N < B and thus (4.29) and (4.30) are true. The only way = 0 is if T < B everywhere in fl. Thus for B > 1 V3 the whole region fl is unyielded and the bubble does not move. 42 (4.35) Comparison with the First Stopping Condition Using our first stopping condition, B Z ^ - J ^ r W n - W Jz_ - V * ) , (4,6) f \ 2 , i ) 2 we can also get an estimate for the critical Bingham number, above which the bubble will not move. For a spherical bubble (or any bubble that is symmetric through a horizontal plane) the surface integral in (4.36) is zero. This is most easily seen by recalling that the integral is antisymmetric in z. Thus for any bubble shape which is the same under the transformation z —> —z, the integral must be zero. This leaves us with B>^=(z+-z-). (4.37) (z — Z-) is just the diameter of the sphere which is 2. Thus from the first stopping condition + we get an estimate for the critical Bingham number of (4.38) B > -L Our new estimate of B > is better and is possibly quite sharp. The estimate we get from the first stopping condition is less accurate since we have made many more approximations to obtain the more general result. 4.3 Third Stopping Condition: for Long Cylindrical Bubbles Here we consider the case of long cylindrical bubbles. For our analysis we assume that the bubble is formed of a long cylindrical body with spherical cap ends (see figure 4.6). Furthermore we divide the fluid domain fi into three regions, fii, fi2, and f^; fi being the horizontal section 3 of fluid where the bubble has a cylindrical body, and fii and regions of fluid above and below ^ 3 . 43 are respectively the remaining Q 3 1 Rc Q 2 \ Figure 4.6: A long (r <C L) cylindrical bubble. b We begin with the stopping condition (4.2): For > B SUP (_ f ^ L , , (4 9) the only possible solution is one with u = 0 over fl. We first note that since our bubble is symmetric through a horizontal plane, the surface integral term in (4.39) is zero. Thus (4.39) becomes f us dfl n B > - J. n 3 h i W dfl . We can make the further simplification that fn ^3 dfl / u dfl n J i{u)dfl 3 J j(u)dfl' n n3 So for In i(u)dfl 3 we will also have no flow. 44 (4.40) We first consider the numerator of expression (4.41): — f u dQ. We note that Q - 3 / u dQ = V U Ja 3 b b = Q ^ 3 + Trr I?j 2 b U, (4.42) b where V is the bubble volume, U is the bubble velocity, r is the radius of the cylindrical portion b b b of the bubble, and L is the length of the cylindrical portion of the bubble (see figure 4.6). Also, from our non-dimensional scaling V — (since V = b and R = \j^V ) b o An —r s b 92 t + nr L b b 47r so (4.43) = —, and for I > rj we get 1. or 3r 2' 6 (4.44) Hence - J To obtain an estimate of f n u dQ~nr LU 2 3 b b + (4.45) 0( 7(u) dQ we approximate the flow as a one-dimensional axisym- metric flow. For a long bubble, the effects at the ends of the cylindrical section where this assumption breaks down should not significantly affect the order analysis. The flow in Q will 3 have velocity profile similar to that in figure 4.7. In this case the only non-zero component of bubble surface column wall unyielded r = n yielded R =R C Figure 4.7: Velocity profile of the flow around the cylindrical section of the bubble. the stress is At the bubble wall, where we have zero tangential stress, the fluid will be 45 unyielded. The fluid will only be yielded in a narrow region of thickness a near the column wall. The mean velocity of the flow is f/avg = n{R -r ' )' <il ^° 2 b 2n -L rRc r I 7(11) dfl ~ / / 7n Jo Jo / -——drd(f>dz JR -a dr c 3 TTT U 2 h h <n^nn j { 2 7 T a R c L ) 2-Kr R LU 2 b ~ c b R ? - n • 2 ( 4 4 6 ) Combining (4.45) and (4.46), we see that for r <C L, the stopping condition (4.41) becomes: b For 1 p 2 „ 2 S ^ ^ T " ' (4-47) the only possible solution is one with u = 0 over f2. We can also write this in terms of L, using (4.44) as 1 B D 2 - 2 ^ 46 M 4_ L - ( 4 - 4 8 ) Chapter 5 Parameter Dependence Using the variational inequality from the rate of strain minimization, (3.44), we can determine the monotonicity of the individual dissipation terms with respect to the physical parameters of the problem. The dimensional version of (3.44), suitably simplified, is \faJn I z I ~ lij{u)lij{v n l -u)d£l + T Y ( \Jn / 7 ( « ) d f i - / 7(u)dfH Jn J > ~Pe9 [ (fa-u )d£l Jn 3 - £ [ (4- + 4-J (i>i - Ui)n ds. Jon \Ri R2J bti (5.1) b Note that the physical parameters fa, fy, Pe, and £ (respectively consistency, yield stress, liquid density, and surface tension) all appear individually in front of separate terms, pe appears in front of the dimensional equivalent of the a(-, •) term, f y appears in front of the dimensional equivalent of the j(-) terms, etc. To determine the monotonicity of the dissipation terms with respect to the physical parameters, we would like to be able to vary one parameter leaving all other parameters (including bubble shape) fixed, and determine the effect on the dissipation terms. For the yield stress and the surface tension, there are the respective dimensionless parameters, B and 8 already contained in (3.44). Thus we can vary B and 8 to determine the influence of varying f y and £. For the consistency and the liquid density, based on (5.1), it suffices to add to a dimensionless parameter, which characterizes a change from a typical value of consistency or of liquid density, in front of the appropriate term in (3.44). For example, to characterize a change in consistency we need only consider a dimensionless parameter pi = 47 where p\ is a typical consistency, and the variational inequality Wa(«, v -u)+ j(v) - j(u) > L[v - u) (5.2) and the functional J(v) = -^—u a(v,v)+j(v)~L(v). n +1 (5.3) i Similarly, to characterize a change in liquid density we use the dimensionless parameter pe = |f, where p\ is a typical liquid density, and the variational inequality a(u,v - u) + j(v) - j(u) > -pi [ (u - u ) dQ - [ 6 (~ + -^-) (^ - Ui)n ds. Jan v-Ki -"-2/ 3 3 bji (5.4) b Generally if the bubble velocity increases so will the dissipation rate terms. Hence from this parameter dependence we are also able to get an idea of how the physical parameters affect the bubble velocity. 5.1 Consistency Consider two consistencies p^ let > pfp. Let be the solution corresponding to and be the solution corresponding to p^p. Now the true solution minimizes the functional (5.3). Also for the true solution 71 J{u) = So for the case for p^\ ^p a{u,u). (5.5) t we have n -p^a( ^,u^) n +1 = -±-p^ (u^,u^) U 1 + j(uW) - L(u^). a n+1 (5.6) (2) Similarly, for the case of p\ we can write a(u<U<>) =n-+^1? W U ) + i ( u ) 2 n -pf n +1 n (22)1 , coi2 n +1 t 2 2 ( 2 ) -/4 a(u( \u ) <. - nL1+- 1^ f a ^ d ) , ^ (2) ( 2 ) 1 -L(uW) ) ) + j(uW) - L(uW). (5.7) Subtracting (5.6) from (5.7) ^^«( >u ) - -^?W \«< >) < n + (1 (1) n (21 2 , to\ 2 COIN U 48 , 1 _^(D^d) j U (D). ( 5 . 8 ) Since ^ > (5.8) implies J^vPa(uM,uW) n + 1 e - -^/i< a(u( U n+1* 2) ' v 2 (2) v ) <0 n+1 n+ 1 c Thus /^(^W^^W^); =» tf^M? 5 (5-9) uga{u, u) is a decreasing function of [ig. 5.2 Yield Stress (Bingham Number) Consider two values of the Bingham Number, B^ ing to B^\ and let be the solution correspond- 2 with v = u^> we get - „(D) + ( [ j(uW) dii - [ ( w ) dii) > L(u& - u^) \Jn Jn ) (1) 7 Using (3.44) for (B^ \ u^), now with v = 2 a(uW, 2 be the solution corresponding to B^ \ Using the variational inequality (3.44) for (B^\u^), a(u^\v.™ < B^ \ Let - we get + B& f [ ( ( D ) dii - ! ( « \Jn Jn 7 (5.10) U 7 ( 2 ) ) dii) > L(u^ J - u&) (5.11) Adding (5.10) and (5.11), and recalling that a(-,-) is linear in its second argument and that L(') is linear a(u^,u^-u^)-a(u^\u^-u^) + (B^-B^)( f -y(u^)dil) > 0. Jn J (5.12) [ j(u^)dil\Jn Using the result from convexity and Gateaux differentiability, (3.51), we can write — n +al( u ( U ( 2 ) —afuf'U n+1 1 1 1 v 2 ) n- —+ «l ( i i ; v ( 1 U ( 1 ) )/ -> va{vP-\uW - ) - -^-a(<z( U ) > a(uW,uM ' n+ 1 / - v 2 v 49 (2) > - ( )) 2 u v(5.13) i (5.14) Adding (5.13) and (5.14) we see that 0 > a(u^\u^ - « « ) - a(u< \ u& - „(D). (5.15) 2 It should be noted that the result in (5.15) is a general result. Combining (5.15) with (5.12), _ 5(2)) ( f ( ( > ) dfl 7 u Thus rffi)> 0. [ (w^)) 2 7 Jn Vin (5.16) ) < 5(2) implies I i{ W)dfl< in f (u in U ( 1 ) 7 )ffi. (5.17) The yield stress dissipation is a decreasing function of the Bingham number. 5.3 Density Consider two liquid densities pf^ < p\ \ Let be the solution corresponding to p^\ and let 2 tt( ) be the solution corresponding to pf \ Using (5.4) for (p^\u^>), with v — 2 a(uM,uW - w ( 1 ) )+i(^ ))-i(^ ))>-pj 2 1 f (4 -4 )df2 2) 1 ) in ian 2 e a(u( ),tt 2 -« (1) ( 2 ) ) + now with v = - J> ( 2 ) 1} 8(i- - f Using (5.4) for (p \u^), we get \m 6 + -U (<> - tiJ^K, d a . 2 (5.18) u M2J we get ) > -Pf / -u?)dfl Jn "Sen P{jk jk) ^ ~ ^ ' + iU U )nb idS ' (5 19) Adding (5.18) and (5.19), and recalling that a(-, •) is linear in its second argument a(uW,uW - u W ) -a{uV\uW - W) > (p[ - p™) [ (u in 2) U {2) - uf) dfl. (5.20) Combining (5.15) with (5.20) we obtain {P[ -P[ ) 2) [ 1] (4 -4 )dn<o. 2) 1} (5.21) 7n Thus with p[ l) < pf\ [ u in {2} dfl < [ u f dfl. in 50 (5.22) This result can be related to the bubble velocity, U , through J^u^dQ = —V U , where V is b b b b the bubble volume. Thus we conclude that the bubble velocity is an increasing function of the liquid density: pU<pM 5.4 => U < Uf\ W b (5.23) Surface Tension Consider two values of surface tension 8^ < 8^ \ Let be the solution corresponding to 2 B^\ and let be the solution corresponding to B^ \ Using the variational inequality (3.44) 2 for (B^\u^), with v = we get a(uW,uW - W)+j( W)-j(uW)>u [ u (u -u^dQ {2) Jn 2J Jan \ i K K b Using (3.44) for {B^ \u^), now with v = we get 2 a(uW,uM -uW)+j(uV)-j(uW)>- [ (4 -u )dQ :) ( 2) 3 Jn -e [ (w {2) +w ) ^ - ^ ^ - - (5 25) Jan b Adding (5.24) and (5.25), and recalling that a(-, •) is linear in its second argument a(uU,uW - (U) - a(uM,uW - (D) > fpW - 0(D) f U (i_ + ±-) U Jdn b \R\ („« - u^)n ds. bil -n-2/ (5.26) Combining (5.15) with (5.26) we obtain (/3 - /3 ) / (2) (1) K + Jan \i K b 2j R (2) - ^K* ds < 0. From our earlier results with the surface tension; see (4.14), and the fact that (5.27) u dQ, = —V U , 3 b b we have that 1 1\ f . -s- + S-) W , i ds = 2nU / f'(f"f -±U. z+ an \Ri b b R-2J J_ z 51 J - (f) - 1) ^ > > dz ((y/)2 2 3 + 1 ) 2 (5.28) Combining (5.28) with (5.27) we get 2n(3^ - P^)(U (2) b ~U ) [ {1) Z+ b dz < 0. ^ ^-^ ) - ) , , 2 1 s ((/0 + i ) (5.29) 2 J*- 5 It turns out that for "typical" bubble shapes /*+ • H / " / - ( / ' ) - ) dz < 0. Thus ((/') +i) 2 2 1 5 ul <U ; L) BW<BW {2) B the bubble velocity is an increasing function of the surface tension. 52 (5.30) Chapter 6 Experiments 6.1 6.1.1 Experimental Setup Experimental Apparatus All the experimental work cited in the introduction has essentially been done using apparatus of similar construction. Our apparatus is designed along similar lines. Our apparatus consists of a 6' long clear acrylic column that has an inner diameter of 1.5" and an outer diameter of 2". The column is attached to a clear acrylic base. Our viscoplastic fluid is pumped in through the sides of the base of the column using a progressive cavity positive displacement pump. This type of pump was chosen at it does not cause high shear rates in the fluid and it does not result in pulsatile flow. These two features are important in order to avoid damaging the long polymer chains in the viscoplastic fluid. It is the cross-linking of these polymer chains that give the fluid its yield stress. The air bubbles are injected from the centre of the base of the column by a pneumatic cylinder. The cylinder is controlled by a linear actuator which in turn is controlled through a computer. By controlling the linear displacement of the cylinder we can control the volume of air being injected into the column. Lastly, a reference grid is located directly behind the column. Each experiment, consisting of the injection of one bubble into the column and observation of its progress after reaching a steady state, is filmed using digital video. Assuming an axisymmetric profile for the bubble, we can extract all the necessary data from the video afterwards (see Appendix C). A schematic of the experimental setup is shown in figure 6.1. 53 clear acrylic column video 0 era) camera pneumatic cylinder pump 1 1 | | controller / mechanical actuator 4 X Figure 6.1: Schematic diagram of the experimental setup. 6.1.2 Viscoplastic Fluid - Carbopol The viscoplastic fluid used in our experiments is Carbopol. It is a rheology modifying polyacrylic based polymer made by Noveon Inc. (formerly B.F. Goodrich Co.) Carbopol is an ideal experimental fluid in that it is clear and comes very close to performing as an ideal viscoplastic fluid, without any elastic behaviour at low shear rates [34]. Two Carbopol polymers (940 and EZ2) are used in order to obtain a wider range of fluid parameters. In general, the Carbopol EZ2 has a lower yield stress than Carbopol 940 for a given concentration. And for both Carbopols the yield stress increases with concentration. Carbopol can be appropriately modelled as a Herschel-Bulkley fluid. The rheological parameters of the Carbopol solutions (yield stress, fy; consistency, jlf, and power law index, n) are measured using a cone and plate rheometer (Bohlin Instruments). The yield stress is directly measured by successive creep/recovery tests. This is the method suggested by Coussot [35]. A constant stress is applied to the fluid and the strain response after a fixed time is measured. When the strain response is plotted versus the applied stress there is a stress value above which the strain 54 response begins to increase significantly as a function of the applied stress. The stress value at which this occurs is defined to be the yield stress. The consistency, fa, and power law index, n, are determined from a controlled stress viscometry test, where the rate of strain is measured at various stresses. The parameters fig and n are fitted to the data in a least-squares sense. For a more detailed description of the parameter measurements and an example calculation refer to Appendix D. The density of the fluid is measured using standard specific gravity bottles and a mass balance. 6.1.3 Preparation of C a r b o p o l Solutions Carbopol comes as a fine polymer powder. The powder must first be dispersed in an aqueous solution and then an alkaline neutralizer must be added to thicken the solution and to produce a yield stress. We begin with 10 L of distilled water and, while mixing at a constant rate of 300-500 rpm with a Heidolph medium-duty laboratory mixer, we slowly add the desired amount of Carbopol. Once all the Carbopol has been added, the solution is mixed for an additional 4-8 hours. This is to allow for the complete hydration and mixing of the powder. The solution is then left to stand for another 4-6 hours to allow all air to escape and to allow any foam to break up. It is important to allow all bubbles to disappear before the solution is neutralized. Once the fluid has a yield stress it becomes exceedingly difficult to eliminate bubbles from the fluid. The fluid is neutralized with a sodium hydroxide solution. While mixing at a lower rate (~ 100 rpm), the sodium hydroxide solution is added about a millilitre at a time. In between the solution is left to mix for several minutes. It is important that the rate of mixing is low enough that air is not entrained into the Carbopol. As the sodium hydroxide is added the consistency of the Carbopol increases and a yield stress develops. Thus the rate of mixing can slowly be increased with the addition of the sodium hydroxide. Sodium hydroxide is added until the Carbopol reaches a pH of 7.0. This requires roughly a 2:1 ratio (by mass) of sodium hydroxide (20% solution) to Carbopol powder. 55 To obtain a higher density solution sugar is first added to the distilled water. 6.1.4 Experimental Method Prior to each experiment the column is drained, the fluid cycled through the pump, and then the column is refilled. This is to ensure that the fluid is adequately mixed between experiments. Gheissary and van den Brule [15] conducted experiments to examine the time effects and wake effects on spherical particle settling in a Carbopol solution. They observed the following: For very short time intervals between spheres begin released in the Carbopol (approximately 2 seconds) there was no noticeable change in the settling velocities. For time intervals on the order of ten minutes the settling velocity increased for following spheres to some higher limiting value. For longer time intervals the limiting settling velocity increased and the number of successive drops required to reach that value decreased. Furthermore, after waiting two hours, in an attempt to let the fluid "recover", the settling velocity of the sphere increased even further. It required approximately 72 hours for the fluid to restore itself to the initial state. Alternatively, mixing the fluid also seemed to restore it to its initial state. This strongly time dependent fluid behaviour requires one to take great care to ensure the fluid is well mixed before the column is filled. Furthermore, we also allow the bubble to travel ~ 50 cm before any measurements are taken. This is to allow the bubble to reach steady state. In the worst case this distance is approximately equal to 10-12 times the effective bubble radius, R, which is still more than adequate for the bubble to reach steady state. There is also another interesting complication that arises because of the yield stress of Carbopol. In the unyielded regions all that we are able to know about the stress is that it is below the yield stress. In reality the stress in the unyielded region is non-unique (the velocities will be unique though); for a given velocity field in fi, we can have different stresses in the unyielded region. For this reason we cannot accurately control the volume of the injected bubble, despite being able to accurately control the cylinder actuation. Since the unyielded fluid can support a range of stresses, the injection pressure can vary over a small range. Thus the volume of the actual bubble will vary in a small range around the injected volume. For small bubbles this 56 could result in as much as 35% variation in volume. To prevent this from causing complications we instead calculate the volume of the bubble from the video, again making the assumption of an axisymmetric bubble profile. This leads to an accurate estimate of the true volume (see Appendix C). 6.2 Experimental Results For our experiments we use six different mixtures of Carbopol. The mixtures and their fluid parameters are listed in Table 6.1. In each fluid we observe a series of bubbles of various volumes Consistency, p,g Yield Stress, fy Carbopol Density, p Power Law Solution (kg • m~ ) Index, n (Pa • s ) (Pa) 1 1068 0.37 4.3 2.3 2 1067 0.35 6.8 4.0 3 1066 0.36 7.5 5.5 4 1158 0.39 5.8 2.5 5 1189 0.46 2.5 5.2 6 1072 0.42 3.2 2.2 t z n Table 6.1: Fluid properties of the Carbopol mixtures. ranging from approximately 0.5 mL to 50 mL, and measure their velocities. The general trend for a given fluid is that as the volume increases the velocity increases rapidly; then with a further increase in volume the bubbles begin to fill the cross-section of the column, and the velocity levels off at some constant value (see figure 6.2, solution 6). In some cases the velocity will increase to some maximum value and then decrease before it levels off to a constant value (see figure 6.2, solution 2). The constant velocities at which sufficiently large bubbles > 0.5^ travel in each of the solutions is given is Table 6.2. The possible variables upon which the bubble velocity could depend are bubble size, yield stress, fluid density, gravity, consistency, column radius, and possibly surface tension and the power 57 18 . r ±_ . 16 14 12 solution 2 solution 6 §10 { V 10 20 30 volume (mL) 40 50 F i g u r e 6.2: Curves of velocity versus volume for two C a r b o p o l solutions. law index (R, f y , pg, g, jig, R , c £, and n respectively). I n c l u d i n g the b u b b l e velocity, Ub, (our dependent variable), a n d a p p l y i n g dimensional analysis we o b t a i n six dimensionless quantities: Ub = Fr, Froude number, IgR ~^=B, pegR Pl9 Pe. B i n g h a m number, = Re*, modified R e y n o l d s number, R_ R r and PegR 2 n. It is reasonable to assume that the surface tension is constant for a l l experiments; therefore, from the B u c k i n g h a m P i theorem we have that for our results " U A R where / is some u n k n o w n function. 58 P * R (6.1) Carbopol Solution Velocity of Bubbles with > 0.5) (cm/s) 1 11.1 2 5.5 3 5.6 4 13.0 5 10.5 6 16.8 Table 6.2: Velocities of bubbles for which > 0.5 in the different Carbopol solutions. The complete data set of all bubbles in all solutions is given in figure 6.3. Again we can see the characteristic increase of velocity with bubble size until the bubbles begin to fill the column (]T ~ ^ ' ^ ) ' * ^ ^ P a w c o m t the velocity levels off to a constant. In figure 6.3 a curve of constant velocity is given by y = —t since Fr oc — VR V R Furthermore, for -§- > 0.5, since we have a constant velocity, our system can be considered to Rc have two regimes; one being < 0.5 where / in (6.1) is dependent on -§-, and the other being Rc Rc > 0.5 where / is independent of R c Rc When comparing the six solutions it may appear that the values of n are relatively close; however, even these differences appear to be significant. To see the effect of changing n on the bubble velocities we first note that for solution 1 and solution 6 the Bingham numbers and the modified Reynolds numbers all lie along the same curve; see figure 6.4. However, from Table 6.2 we can see that the velocities for the bubbles with -§- > 0.5 are 11.1 cm/s and 16.8 cm/s for Rc solutions 1 and 6 respectively. This difference in velocity must be the result of the different values of n for the two solutions (0.37 and 0.42 respectively). To examine the dependence of the velocity on the Bingham number and the modified Reynolds number we consider the data from the first four solutions only. Since the values of n for these 59 0.7 0.6 II I 0.5 I fc, a | 3 I solution 1 0.4 solution 2 A solution 3 x solution 4 i § 0.3 • solution 5 o solution 6 a• Hp 0.2 • o™ . • •• i 0.0 0.2 0.4 0.6 0.8 Bubble Length 1.2 1.4 1.6 # Figure 6.3: The complete data set. Froude number (non-dimensional velocity) plotted versus non-dimensional bubble length. solutions differ by only about 10%, we hypothesize that the variation due to the differences in n should be small compared to the variation due to differences in the Bingham number and the modified Reynolds number. From the data for the first four solutions we fit, in a least-squares sense, a two-dimensional surface to the data: = f(B,Re*). (6.2) Contours of this surface are shown in figure 6.5. Using this surface we are able to plot crosssections of the surface to examine the dependence of the velocity on the Bingham number and on the modified Reynolds number separately. Some cross-sections with the modified Reynolds number held fixed are given in figure 6.6. And cross-sections with the Bingham number held fixed are given in figure 6.7. From figure 6.6 we see that for fixed Re* as we increase the Bingham number the velocity 60 25 " o solution 1 solution 6 •o20 §15 ;g o c |lO 0.02 0.04 Bingham Number, B 0.06 Figure 6.4: Modified Reynolds number plotted versus Bingham number. For the two solutions the data lie along the same curve (dotted line). decreases. This makes sense physically in that if all other variables are held fixed and we increase the yield stress (increasing the Bingham number) then the bubble velocity should decrease. The results from figure 6.7 are more surprising. We see that as the Re* decreases from a large value the velocity begins to decrease; however, as Re* continues to decrease the velocity reaches a minimum and then begins to increase rapidly. Normally we would expect that if all other variables were held fixed and we increased the consistency of the fluid (decreasing Re*) the drag on the bubble would increase causing the velocity to decrease. For the larger modified Reynolds numbers this appears to hold true, however at lower values this is no longer the case. It is not clear what physical phenomena produce this result. 6.3 6.3.1 Comparison with Analytic Results Parameter Dependence With regards to the parameter dependence, the experiments verify that an increase in the yield stress (or Bingham number) results in a decrease of velocity (result (5.17) from section 5.2). Our result for the consistency (5.9), however, only appears to be valid over certain regions of values for fig. Also for an increase in density, which results in B decreasing and Re* increasing, 61 0.01 0.02 0.03 0.04 0.05 0.06 0.07 Bingham Number, B Figure 6.5: Contours of the Froude number, 0.08 0.09 X^, plotted as a function of the Bingham number and the modified Reynolds number. The contour spacing is 0.025. over the range of modified Reynolds numbers where the velocity is an increasing function of Re*, we have an increase in velocity (result (5.23)). Finally since we did not specifically alter the surface tension of the solutions we are unable to verify the result concerning changes in the surface tension (5.30). 6.3.2 Stopping Conditions While it would have been ideal to have data on bubbles that moved as well as bubbles that were stopped, with the current experimental setup it was not possible to examine stopped bubbles. Firstly, it is extremely difficult to inject a stopped bubble into the column that is not still attached to the injection nozzle. And secondly, a stopped bubble is, in general, not axisymmetric and thus it is impossible to calculate the volume from the video data. We must then suffice with only having data on moving bubbles with which to compare our stopping condition results. 62 It should also be noted that the stopping conditions are independent of the consistency and power law index, Re* and n. While they influence the flow for moving bubbles, if a bubble is stopped in a given fluid, changing the consistency or power law index will not suddenly cause the bubble to start moving. The stopping conditions are only functions of the geometry of the problem (bubble shape and column radius). (Though, indirectly through the nondimensionalization, they depend on the density and gravity.) Our first stopping condition (4.26) is B > ±. - )(l- ^lf 2 (z+ - 2 ^ ~^ Z+ z + ' y v b j_ z -V d z ((fy+ To evaluate this it is necessary to know the value for the surface tension between Carbopol and air. Unfortunately for a yield stress fluid it is not easy to measure the surface tension without effects of the yield stress interfering. Conventional methods using capillary tubes or the Du Nouy ring method will not be able to separate surface tension effects from yield stress effects. (It should be noted that Kim et al. have devised a method of measuring the surface tension of a yield stress fluid [36].) Fortunately for our analysis, a precise value of the surface tension, £, is not necessary. Using an approximate value of £ = 0.0725 N/m (the value of the surface tension between water and air), the value of the surface integral term in (6.3) turns out to be 63 1 1 0.4 -- B B B B = = = = 0.02 0.03 0.04 0.06 co O Fro ude Number, LL \ \ \ \ \ \ \ \ \ \ \ \ \ i V - - . \ / 0.1 I 0 I I I 5 10 15 I 20 Modified Reynolds Number, Re Figure 6.7: Froude number as a function of modified Reynolds number for fixed B. several orders of magnitude smaller than the other term, i.e., 27T/3 /b J /•*+ / ' ( / " / - (/') r -1) 2 ni" i - ( / y - D , D * Jz- < < L ( M ) ((/') (( f'\2+ l ) 2 In fact for our data 2*0 \r* ^rwi -0.0008 < —fV J_ ( L b L - ^—L—_—Ldz< in -1) 0.001. 2 1 z ( / 0 2 + 1 (6.5) ) 2 Assuming the true value of £ to be of the same order as that for water and air, we simplify (6.3) for our following analysis: B ~ 2^ ~ - ' {Z+ Z ] In figure 6.8 we plot our data with non-dimensional bubble length, z ( 6 + ' 6 ) — z_, versus Bingham number. The solid line represents the curve B — 7^{z+ — z-). We can see that all the bubbles (none of which were stopped) respect condition (6.6). While none of the bubbles were spherical, the smaller bubbles could, as a very rough first approximation, be considered to be spheres. We can then use our second stopping condition (4.35) to see whether it is still a valid stopping condition. The dotted line in figure 6.8 represents the curve B = Again we can see that all the bubbles respect condition (4.35): B > ± 64 (6.7) Figure 6.8: Experimental results plotted with non-dimensional bubble length versus the Bingham number. The solid line represents the curve B = ^^i + ~ -) d the dotted line represents the curve B = -4= z z a n While the bubbles with (z — zJ) small are closer to spheres in shape it appears that condition + (6.7) is applicable to all the bubbles. Though condition (6.7) might not hold in general, one reason it may appear this way is that as (z — zJ) becomes larger the bubbles tend to move + with a higher velocity than the smaller bubbles. Thus these bubbles are "farther away" from stopping and hence condition (6.7) could appear only to be valid. Also, if we plot the data with the Froude number versus the Bingham number, and compare this to condition (6.7) (see figure 6.9), we see that for increasing Froude numbers the bubbles are farther away from the line B = i.e., are "farther away" from stopping. Finally we have our third stopping condition for long cylindrical bubbles, (4.47): B - 2 ^ i ^ - ( 6 - 8 ) Since none of our bubbles are actually cylindrical they do not have a constant value of r . Nevb ertheless we would like to be able to roughly compare this condition with our other conditions and with our data. So with that interest in mind we make the approximation r 2 b 65 = ~ ? — ^ — r ; see 0.6, X 0.5 x X x i t 0.4 X X x cu x JZI E X* >\* • CD - X x x %* X x X x x • x x x*<* x *x X * * 0.1 - X x X X X x Xx 0 * 1 Bingham Number, B Figure 6.9: Experimental results plotted with Froude number versus the Bingham number. The dotted line represents B = Note that B decreases with increasing Fr, and thus faster bubbles are further from the B = -4= line. (4.44). This represents a sort of average effective radius of the bubble, and thus (6.8) becomes 1 ^ - 3 ( ^ 7 ~ 2 R V c ' Since the data is scaled according to the length scale R, which is different for each bubble, it is difficult to plot condition (6.9) as a function of (z — zJ). Instead we separately calculate the + critical Bingham number for each bubble and plot this along with the data (see figure 6.10). The dashed line in figure 6.10 is a best-fit line through the critical Bingham number data. Again all the data appears to respect condition (6.9). Unlike our other two conditions, as the bubble length increases the critical Bingham number decreases. This represents the idea that if a bubble is very large it can more easily fill the column entirely such that there is no layer of fluid between the bubble wall and the column walls. In this case the bubble will no longer be able to move; it has essentially separated the fluid region into two disjoint regions where fluid cannot pass from one region to the other. Using this best-fit line we can compare all our stopping conditions to the data. Figure 6.11 shows all the data plotted with the non-dimensional bubble length versus the Bingham number 66 —i \ * 1 1—'»«•'! x . \ X 0) XX * \ X x # - " - •\ •. x " x"x \ \ \• \ .\*. • •\- .r V' . \-> . X X *x X * ." - i\ **. *\ , X N - Ql ' 10 10 1 -2 B cnt • \ \ • data V' •< . -1 10 Bingham Number, B 10 0 ' 1 10 Figure 6.10: Experimental results plotted with non-dimensional bubble length versus the Bingham number. The critical Bingham numbers for the data are also plotted (dots). The dashed line is a best fit line through the critical Bingham numbers. with all three stopping conditions included. 67 10 10 Bingham Number, B 10 Figure 6.11: Experimental data combined with all three stopping conditions. 68 Bibliography [1] J.G. Oldroyd. A rational formulation of the equations of motion of a plastic flow for a Bingham solid. Proceedings of the Cambridge Philosophical Society, 43:100-105, 1947. [2] O.L.A. Santos and J.J. Azar. A study of gas migration in stagnant non-newtonian fluids. Technical report, Society of Petroleum Engineers, 1997. SPE 39019. [3] A. Johnson, I. Rezmer-Cooper, T. Bailey, and D. McCann. Gas migration: Fast, slow or stopped. Technical report, Society of Petroleum Engineers, 1995. SPE/IADC 29342. [4] A.B. Johnson and D.B. White. Gas rise velocities during kicks. Technical report, Society of Petroleum Engineers, 1990. SPE 20431. [5] S.M. Bhavaraju, R.A. Mashelkar, and H.W. Blanch. Bubble motion and mass transfer in non-Newtonian fluids: Part 1. single bubble in power law and Bingham fluids. AIChE Journal, 24(6): 1063-1070, 1978. [6] S. Stein and H. Buggisch. Rise of pulsating bubbles in fluids with a yield stress. Z. Angrew. Math. Mech., 80(ll-12):827-834, 2000. [7] A.N. Beris, J.A. Tsamopoulos, R.C. Armstrong, and R.A. Brown. Creeping motion of a sphere through a Bingham plastic. J. Fluid Mech., 158:219-244, 1985. [8] D.D. Atapattu, R.P. Chhabra, and P.H.T. Uhlherr. Creeping sphere motion in HerschelBulkley fluids: flow field and drag. J. Non-Newtonian Fluid Mech., 59:245-265, 1995. [9] J. Blackery and E. Mitsoulis. Creeping motion of a sphere in tubes filled with a Bingham plastic material. J. Non-Newtonian Fluid Mech., 70:59-77, 1997. [10] V. Dolejs, P. Dolecek, and B. Siska. Drag and fall velocity of a spherical particle in generalized newtonian and viscoplastic fluids. Chemical Engineering and Processing, 37:189-195, 1998. [11] G Saha, N.K. Purohit, and A.K. Mitra. Spherical particle terminal settling velocity and drag in Bingham liquids. International Journal of Mineral Processing, 36:273-281, 1992. [12] B.S. Padmavathi, T. Amaranath, and S.D. Nigam. Stokes flow past a sphere with mixed slip-stick boundary conditions. Fuild Dynamics Research, 11:229-234, 1993. [13] L. Jossic and A. Magnin. Drag and stability of objects in a yield stress fluid. AIChE Journal, 47(12):2666-2672, 2001. [14] J. L i and Y . Y . Renardy. Shear-induced rupturing of a viscous drop in a Bingham liquid. J. Non-Newtonian Fluid Mech., 95:235-251, 2000. [15] G. Gheissary and B.H.A.A. van den Brule. Unexpected phenomena observed in particle settling in non-Newtonian media. J. Non-Newtonian Fluid Mech., 67:1-18, 1996. 69 R. Clift, M.E. Weber, and J.R. Grace. Bubbles, drops, and particles. Academic Press, 1978. Z. Zapryanov and S. Tabakova. Dynamics of bubbles, drops and rigid particles. Fluid Mechanics and its Applications. Kluwer Academic Publishers, 1999. F.P. Bretherton. The motion of long bubbles in tubes. J. Fluid Mech., 10:166-188, 1961. S.V. Vasil'chenko and A . G . Potapov. Gas bubble dynamics in a vicsoelastic-plastic medium. Heat Transfer Research, 27(l):4-8, 1996. W. Prager. On slow visco-plasticflow,pages 208-216. Academic Press Inc., 1954. Presented to Richard von Mises by Friends, Colleagues, and Pupils. P.P. Mosolov and V.P. Miasnikov. Variational methods in the theory of the fluidity of a viscous-plastic medium (Variatsionnye metody v teorii techenii viazko-plasticheskoi sredy). PPM, 29(3):468-492, 1965. P.P. Mosolov and V.P. Miasnikov. On stagnant flow regions of a viscous-plastic medium in pipes (O zastoinykh zonakh techeniia viazko-plasticheskoi sredy v trubakh). PPM, 30(4):705-717, 1966. P.P. Mosolov and V.P. Miasnikov. On qualitative singularities of the flow of a viscoplastic medium in pipes. PPM, 31(3):581—585, 1967. R. Glowinski. On the approximation of an elliptic variational inequality of the Bingham type. Revue Francaise d'Automatique Informatique Recherche Operationnelle, 100:13-30, 1979. R. Glowinski. Numerical methods for nonlinear variational problems. Springer series in computational physics. Springer-Verlag, 1984. R.R. Huilgol. Variational principles and variational inequalities for a yield stress fluid in the presence of slip. J. Non-Newtonian Fluid Mech., 75:231-251, 1998. R.R. Huilgol and Q.D. Nguyen. Variational principles and variational inequalities for the unsteady flows of a yield stress fluid. International Journal of Non-Linear Mechanics, 36:49-67, 2001. R.R. Huilgol. Variational inequalities in the flow of yield stress fluids including inertia: Theory and applications. Physics of Fluids, 14(3):1269-1283, 2002. LA. Frigaard, O. Scherzer, and G. Sona. Uniqeness and non-uniqueness in the steady displacement of two visco-plastic fluids. Z. Angrew. Math. Mech., 81(2):99—118, 2001. LA. Frigaard, S. Leimgruber, and O Scherzer. Variational methods and maximal residual wall layers. J. Fluid Mech., 483:37-65, 2003. M.P. do Carmo. Differential geometry of curves and surfaces. Prentice Hall, Inc., 1976. H. Lamb. Hydrodynamics. Dover Publications, sixth edition, 1945. 70 [33] G.K. Batchelor. An introduction to fluid dynamics. Cambridge University Press, 1967. [34] S.J. Curran, R.E. Hayes, A. Afacan, M.C. Williams, and P.A. Tanguy. Properties of Carbopol solutions as models for yield-stress fluids. Journal of Food Science, 67(1):176180, 2002. [35] P. Coussot. Mudflow rheology and dynamics. IAHR monograph series. A.A. Balkema, 1997. [36] S. Kim, Y.I. Cho, W.N. Hogenauer, and K.R. Kensey. A method of isolating surface tension and yield stress effects in a U-shaped scanning capillary-tube viscometer using a Casson model. J. Non-Newtonian Fluid Mech., 103:205-219, 2002. 71 Appendix A A Result on the Effect of Walls If we consider the flow of a bubble in an infinite domain, there will be a finite distance away from the bubble where the fluid will be unyielded. In this unyielded region far away from the bubble, we can "disturb" the fluid without it affecting the propagation of the bubble. So for example, we can introduce walls (of any shape) anywhere in this unyielded region and the flow in the region surrounding the bubble will not change. Or if there are already walls in this unyielded region we can move the the walls within this unyielded region without affecting the flow. Physically this makes sense in that to calculate the flow around the bubble the boundary conditions (not on the bubble surface) must be applied at the yield surface. Hence if we disturb the fluid outside the yield surface without it affecting the shape of the yield surface, there is no way for that information to reach the bubble. Mathematically this results follows from our rate of strain minimization (3.53): The true velocity field minimizes the functional ' Ja(v) = -^-a(v,v)+j(v)-L(v) (A.l) n +1 over the functional space V(tt), where V{tt) = | u = («i,« ,«3) 2 Vi = 0 i n fi, a n d ^ = 0 on dfi^j e C°°(fi), (A.2) and all the integrals in (A.l) are over tt. Consider two domains fii and $7.2 with tt\ Ctt,2-In each domain we have a bubble of exactly the same shape, that is dfii^ = dtt.2^- We know that if a velocity field minimizes (A.l) then it is the unique minimizer. So we let minimize Jn>i(w) over V(fii) and let minimize Jct ( ) V(tt,2). Furthermore we assume that in fi2 \ {fii} = 0, i.e., that the walls of both domains are in the unyielded region. This arrangement is depicted infigureA . l . v o v e r 2 On dtti (the walls), = 0 and on dtti^, satisfies the necessary boundary conditions, namely (3.20)-(3.22). Thus, truncated appropriately, G V(tti). tW Since in tt \ {tt }, We claim now that 2 x =0 <y( W) = 0 in tt \ {fii} => J n ( u ) = 0 in tt \ {fii}. also minimizes Jcii(v) over V(tti). If it did not there would exist (2) U 2 72 2 2 unyielded region fi 2 0,1 2,6 2,ixi / Figure A . l : A bubble propagating in two different domains. u* e V(Qi) such that Jn^tt*) . < f and the velocity field Jni(u^) u* in Qi 0 in n \ {fii} ' 2 would give Jn (<0 < 2 Jn (u 2 ( 2 ) ) which is impossible. Therefore minimizes Jo^ («) over V ( f i i ) and by uniqueness in fii. The flow does not depend on the position of the walls in the unyielded region. 73 Appendix B Some Differential Geometry Results for Surfaces These results are taken from do Carmo [31]. In general, for a two dimensional surface parameterized by x = x(s,t), where Si < s < sj and U <t < tf, the metric tensor, g p, for the surface is given by a *»(•.«)=&•-: where x = s *;.:;)• and x = f j r . t A unit normal to the surface is given by N(s,t) = X t X X (B.2) s x\ \Xf X s And the second fundamental form tensor is given by where x ss = f^f, x st = | ^ | , and x = ^f. tt Then the mean curvature, H(s,t), of the surface is given by And for an integral over the surface / ds = / Jx(s,t) J Si / ||x x x \\dtds. s t (B.5) J ti Now for a "bubble-like" surface of revolution, revolved about the z-axis, as in figure B . l , one possible parameterization is given by x(z,6) = {xi,x ,x ) 2 3 = (x,y,z) = •(/(*) cos 0,/(z) sin0, z), 74 (B.6) Figure B . l : Surface of revolution. where z_ < z < z and 0 < 0 < 2ix. So + cc = (/'cos0,/'sin0,l), z ^ 0 = (-/sin0,/cos0,O) and thus the metric tensor, (B.l), is given by (z,6) 9aP =[ { f ) l + (B.7) £ 1 Note that there is no 0 dependence in the metric. The tensor g @ is simply the inverse of g p: a a 0 U y (B.S) 2 From (B.2) we obtain the outward unit normal for the surface: 1 N = fV(f') 2 (/ cos 0, / sin 0, - / / ' ) . (B.9) +1 From our parameterization (B.6) x = (/"cos0,/"sin0,O), XzO = (-/'sin0,/'cos0,O), = (-/cos0, -/sin0,O), zz and therefore the second fundamental form tensor, given by (B.3), is o 75 (B.10) T h u s the m e a n curvature of our surface of r e v o l u t i o n is given b y R*J 2 V-Ri 2 \((/') 2 +i) fV(f') 2 f +V ' and we have the result i n Ri / " / - (/') - 1 2 Ri) /((/')2 + i ) 2 T h e sign of the mean curvature depends o n the choice of the n o r m a l . If we h a d chosen the i n w a r d n o r m a l of the bubble, there w o u l d be a negative sign i n front of the right h a n d side of ( B . l l ) . T o determine w h i c h sign we need, consider the case of a c y l i n d r i c a l b u b b l e where f(z) — C , a positive constant. T h e n from ( B . l l ) we o b t a i n Ri R2) C + However we want the curvature of the c y l i n d r i c a l b u b b l e to be positive. Therefore we must a d d a negative sign to the right h a n d side of ( B . l l ) . T h u s we have the final result ^ Ri R2J + ,_/•/-(/•)•-1, /((/o + i) 2 ( B . 1 2 ) 2 A l s o , for a n integral over this surface of revolution, we have p2-7r f ds= Jx{z,6) [* JO r+ z [ f^{f') JzZ+ 2 76 + ldzd9. (B.13) Appendix C Data Extraction and Error Analysis C l Velocity Calculation From the digital video of each experiment we calculate the velocity by measuring the time required for a bubble to travel a fixed distance. The video is taken at a rate of 60 frames per second. Figure C l shows two frames from the video of an experiment. Using the reference Figure C l : Video frames of an experiment, (a) t = Os, (b) t = 1.28 s grid (spacing of 0.2") we determine that the bubble has travelled a distance of 8.79 cm ± 0.5%. This represents a measurement error of ±0.1 squares. From the frame indices on the video we determine the elapsed time between the two frames to be 1.28 s ± 0.7%. This represents an measurement accuracy up to ± ^ Q S. Thus, in this case, the velocity is 6.87 cm/s ±0.9%. The 77 error calculation is done separately for each experiment, as the percentage error will vary with the distance travelled and time between the given frames. It should also be noted that the scaling effect of having the reference grid behind the bubble is accounted for when calculating the distance the bubble travels. C.2 Shape Dependent Quantities While the velocity can simply and easily be calculated by examining two frames and determining the distance the bubble has travelled in a specified amount of time, for quantities that depend on the shape of the bubble (volume, surface integrals) we need to extract the shape of the profile of the bubble from the video. In order to do this we apply an edge detection algorithm (Matlab's "edge" function with the Laplacian of Gaussian method) to an image of the bubble. With a little cleaning up of the result to eliminate false edges, we can obtain the profile of the bubble. Figure C.2 shows the actual bubble and the profile obtained from the edge detection. Again the effects of scaling and optical distortion due to the cylindrical geometry (see appendix E) are accounted for in this process. (a) (b) Figure C.2: Profile of a bubble, (a) actual bubble (same as figure C.la), (b) the profile obtained using our edge detection method. The apparent discrepancy in the profile when compared to the video is a result of the cylindrical geometry, for which we have corrected. Once the bubble profile is obtained the volume can be easily calculated; however calculation of the surface integral is still a bit tricky. Due to the finite pixel size the bubble profile is not very smooth when we consider its derivatives with respect to z. Since the surface integral is dependent on the first and second derivatives of the profile with respect to z, we need to 78 somehow smooth the bubble profile to obtain a reasonable result. We do this by fitting a polynomial, in the least-squares sense, to the data for the profile (using Matlab's "polyfit" function). The degree of the polynomial is chosen to be between 6-9. We find that this provides a smooth approximation to the bubble profile without being so accurate as to contain the small oscillations which are the result of the pixelation. C.2.1 A s s u m p t i o n of A x i s y m m e t r y To obtain data pertaining to the shape of the bubble we have made the assumption of an axisymmetric bubble profile. To verify that this assumption is reasonable we compare the profiles obtained by looking at the left half and the right half of a bubble. The profiles are plotted in figure C.3. 5ET f(z) (cm) Figure C.3: Profiles of a bubble obtained by looking at the left half (dashed line) and at the right half (solid line). The assumption is quite good except at the nose and tail of the bubble where the edge detection has difficulty resolving the edges. Nevertheless, this results in an error of only 0.07 mL. This represents a 0.5% error in the volume. C.2.2 A c c u r a c y of the Edge Detection M e t h o d To estimate the accuracy of our edge detection method several circles of know diameter were photographed and the edge detection method was used to obtain a profile for each circle. Using these profiles we then calculated the volume of the equivalent sphere (obtained by revolving the 79 circles) and the value of the bubble surface integral J i - ((/') + l ) 2 ! The results of this are presented Table (C.l). For the circles the true value of the S is zero. Circle Diameter (cm) 2 4 8 Equivalent Volume (mL) 4.19 32.51 268.08 Calculated Volume [mL) 4.27 33.66 268.78 Percentage Error 1.9 3.5 0.3 Actual Surface Integral Value (cm) 0 0 0 Calculated Surface Integral Value (cm) -0.0023 -0.0014 -0.0141 Table C . l : Error in the calculation of shape dependent quantities for test circles. 80 Appendix D Measurement of Rheological Parameters All the measurements are taken using Bohlin Instuments' CVOR 200 rotational rheometer. The tests are done at a constant temperature of 22 °C, which is approximately the temperature of the laboratory in which the experiments are conducted. D.l Yield Stress To determine the yield stress, f y , we conduct a series of creep tests on a sample of the viscoplastic fluid. For each creep test a constant stress is applied to the sample and after 300 s the strain is measured. We normally increase the applied stress between successive tests by 0.5 or lPa. A graph showing the strain after 300 s for various applied stresses is given in figure D . l . We can clearly see two regions, one where the strain does not change significantly when the applied stress is increased, and one where the strain begin to increase significantly with the applied stress. The stress value at which the change between these two regimes occurs (in this case 5.2 Pa) is defined to be the yield stress of the fluid. The measurement of the yield stress can be made to an accuracy of ±0.2 Pa. D.2 Consistency and Power Law Index To determine the consistency, p,£, and the power law index, n , we conduct a controlled stress viscometry test. The rate of strain, 7, is measured at various stresses, T£. This results in the curve given infigureD.2. The constitutive equation for a Herschel-Bulkley fluid can be written as 7(w) = 0 fi{u) if Te(u) < = iJ.a{ii) + f y n fy, if +e(u) > f y . Having already calculated f y , we fit the parameters fit and n to the data in figure D.2 in a least-squares sense. For example, for the sample shown in figure D.2 we obtain the values fae = 2.48 Pa • s and n = 0.462. The masurement of p,g and n can be made to an accuracy of ±6% and ±15% respectively. n 81 appled stress (Pa) Figure D . l : Strain response after 300 s for a sample of Carbopol at various stresses. The lines are to aid visualization. 0 100 200 300 400 rate of strain (s~ ) 500 600 700 1 Figure D.2: Stress plotted versus rate of strain for a sample of Carbopol. The solid line represents the least-squares fit. 82 Appendix E Optical Distortion Due to Cylindrical Geometry An exaggerated top view of the cylindrical geometry of the bubble column is given infigureE . l . r, Ri, and R are respectively the bubble radius (at a given height), the inner radius of the cylinder, and the outer radius of the cylinder. L is the distance from the camera to the front edge of the column and r* is the apparent radius of the bubble, n i , n , and 713 are the indices of refraction of the air, the clear acrylic of the column, and the viscoplastic liquid respectively. Given r*, Ri, R , L, m , n , and n , we need to calculate r. (The other angles, lengths, and 2 2 2 2 3 Figure E . l : Schematic of the cylindrical geometry. points in figure E . l are used to simplify our calculations.) Using the notation of figure E . l , Snell's law applied to the interface between the acrylic column and the viscoplastic liquid is 7J3 sinip = n sintpi. 2 Also from ACDO and ABCO (E.l) 2 respectively, we have r = Ri sin ip (E.2) 2 A= (0 - 0i) + 2 83 fa. (E.3) Combining (E.1)-(E.3) we obtain + 0 ) r = — s i n ( ^ -9i 2 [(sin# cos 6\ — sin#! cos 62) cos 0 2 + (cos# cos#i + sin# sin#i) sin = 2 2 2 fc]. (E.4) n If we can calculate the values of the sine and cosine functions in (E.4), in terms of the given quantities, then we can determine r. 3 Applying Snell's law to the interface between the air and the clear acrylic we have ri2 sin 0 2 From AABO n\ sin <j)\. = (E.5) we see that cj = 0 )1 1 + a . (E.6) Combining (E.5) and (E.6) we have sin 0 2 = — sin(#i + a) n2 Til =—(sin#! cos a + sin a cos #1), Now from 02 = cos AAGO y/l - s i n and (E.7) ^2- 2 (E.8) r tan a = L + R 2 =>sma=— A/T** , + (L + R2) and cosa = —. (E.9) Vr* + (L + R2) 2 2 2 2 In order to calculate sin#i, cos#i, sin# , and cos# from the given quantities, we note that 2 • a 2 sin 0i = —— , and x R2 • a sm9 2 X 2 l A = — Q v R2 ~ x a \/#l cosfli = and 2 7 R2 - 2 cos# = - 2 R\ , , . (E.10) ,. (E-11) 2 2 xi 2 Ri Now all that remains is to determine the values of x\ and x . 2 From AAGO and AABF we can see that r* tana = - L + —, and (E.12) i t .2 tana = (E.13) % 2 t L + R- VR2 - x 2 2 2 2 R — X2 is the distance from the front point of the column to point F. Equating (E.12) and (E.13) and simplifying we obtain a quadratic equation for x : 2 2 2 [{L + R ) 2 2 + r* ] x 2 2 2 - 2r*(L + R ) 2 2 84 x + r* (L 2 2 2 + 2LR ) = 0, 2 which has roots 2r*(L + R ) ± ^/4:r* (L + R ) 2 2 2 4 2 = X 2 - Ar* {L 2 2{(L + R ) + 2LR )[{L + R ) 2 + r* } 2 2 2 2 + r* } 2 2 2 (E.14) Using the requirement that x < r*, we take the smaller of the two roots (i.e., negative square root). Thus using (E.14) with (E.10) we can calculate the values of sin0i and cos0i, in term of the given quantities. Furthermore using (E.9) and (E.10) in (E.7) and (E.8) we can calculate the values of sin <t> and cos fa. 2 2 By considering AFBO and AECO x = x x 2 [\JR we note that for x\ - x 2 2 2 2 - \jRi 2 - xi ) 2 t a n ( 0 i - (f> ). 2 Simplifying we get a quadratic for x\ in terms of known quantities: -2x + 2^R [1 + t a n ( 0 i - <f> )]xi + 2 2 2 2 + x - yR 2 - a ; 2 t a n ( 0 i - fa) 2 2 2 2 - x 2 2 tan(0i - fa) x\ i?i tan (0i - fa) } = 0. (E.15) 2 2 In this case we want the larger of the two roots of (E.15). Note that x is known from (E.14) and tan (0i — fa) can be calculated from the values of sin0i, cos0i, sinfa,and cos fa, which we have already calculated. Thus we can now calculate the value of r. 2 2 85
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Bubble propagation through viscoplastic fluids
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Bubble propagation through viscoplastic fluids Dubash, Neville 2003
pdf
Page Metadata
Item Metadata
Title | Bubble propagation through viscoplastic fluids |
Creator |
Dubash, Neville |
Date Issued | 2003 |
Description | In this thesis we consider the propagation of an air bubble in a cylindrical column filled with a viscoplastic fluid. Because of the yield stress of the fluid, it is possible that a bubble will remain trapped in the fluid indefinitely. We restrict our focus to the case of slow moving or near-stopped bubbles. Using the Herschel-Bulkley constitutive equation to model our viscoplastic fluid, we develop a general variational inequality for our problem. This inequality leads to a stress minimization principle for the solution velocity field. We are also able to prove a stress maximization principle for the solution stress field. Using these two principles we develop three stopping conditions. For a given bubble we can calculate, from our stopping conditions, a critical Bingham number above which the bubble will not move. The first stopping condition is applicable to arbitrary axisymmetric bubbles. It is strongly dependent on the bubble length as well as the general shape of the bubble. The second stopping condition allows us to use existing solutions of simpler problems to calculate additional stopping conditions. We illustrate this second stopping condition using the example of a spherical bubble. The third stopping condition applies to long cylindrical bubbles and is dependent on the radius of the bubble. In addition to our stopping conditions, we determine how the physical parameters of the problem affect the rise velocity of the bubble. We also conduct a set of experiments using a series of six different Carbopol solutions. From the experiments we examine the dependence of the bubble propagation velocity on the fluid parameters and compare this to our analytic results. We find that there is an interesting discrepancy for low modified Reynolds number flows wherein the bubble velocity increases with a decrease in the modified Reynolds numbers. We also compare our three stopping conditions with the data. It appears that all the stopping conditions seem to be valid for the range of bubbles examined despite the fact that when applying the second and third stopping conditions most bubble shapes are not well approximated by a sphere or a cylinder. |
Extent | 4587123 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-10-29 |
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. |
IsShownAt | 10.14288/1.0080989 |
URI | http://hdl.handle.net/2429/14384 |
Degree |
Master of Applied Science - MASc |
Program |
Mechanical Engineering |
Affiliation |
Applied Science, Faculty of Mechanical Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2003-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2003-0400.pdf [ 4.37MB ]
- Metadata
- JSON: 831-1.0080989.json
- JSON-LD: 831-1.0080989-ld.json
- RDF/XML (Pretty): 831-1.0080989-rdf.xml
- RDF/JSON: 831-1.0080989-rdf.json
- Turtle: 831-1.0080989-turtle.txt
- N-Triples: 831-1.0080989-rdf-ntriples.txt
- Original Record: 831-1.0080989-source.json
- Full Text
- 831-1.0080989-fulltext.txt
- Citation
- 831-1.0080989.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-0080989/manifest