UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Macro-size drop encapsulation Maleki, Amir 2014

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

Item Metadata


24-ubc-2014-september-maleki-amir.pdf [ 10.86MB ]
JSON: 24-1.0165999.json
JSON-LD: 24-1.0165999-ld.json
RDF/XML (Pretty): 24-1.0165999-rdf.xml
RDF/JSON: 24-1.0165999-rdf.json
Turtle: 24-1.0165999-turtle.txt
N-Triples: 24-1.0165999-rdf-ntriples.txt
Original Record: 24-1.0165999-source.json
Full Text

Full Text

Macro-size drop encapsulationbyAmir MalekiB.Sc., Sharif University of Technology, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mechanical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2014c© Amir Maleki 2014AbstractViscoplastic fluids do not flow unless they are sufficiently stressed. Whilein some flows this leads to unwanted features, this property can also beexploited in order to produce novel flow features. One example of such flowsare visco-plastically lubricated (VPL) flows, in which a viscoplastic fluid isused to stabilize the interface in a multi-layer flow, far beyond what mightbe expected for a typical viscous-viscous interface. Here we extend thisidea by considering the encapsulation of droplets within a viscoplastic fluid,for the purpose of transportation, e.g. in pipelines. The main advantageof this method, compared to others that involve capillary forces is thatsignificantly larger droplets may be stably encapsulated, governed by thelength scale of the flow and yield stress of the encapsulating fluid. We explorethis setup both analytically and computationally. We show that sufficientlysmall droplets are held in the unyielded plug of the Poiseuille flow. As thelength or radius of the droplets increase the carrier fluid eventually yields,potentially breaking the encapsulation. We study this process of breakingand give estimates for the limiting size of droplets that can be encapsulated.iiPrefaceThis dissertation is an original intellectual product of the author, AmirMaleki with the guidance of his supervisors, Drs. Ian Frigaard and SarahHormozi. The CFD code used in chapters 3-5, originally developed by twoPhD students Andreas Putz and Ali Rostaei, is modified by the authoraccording to the requirement of this study.The results of this study are presented at the “Canadian Applied and In-dustrial Mathematics Society 2014 Annual Meeting” and the “Sixth PacificRim Conference on Rheology”.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . viiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Viscoplastic fluids . . . . . . . . . . . . . . . . . . . . . . . . 21.1.1 Poiseuille flow of a viscoplastic fluid . . . . . . . . . . 51.1.2 Mathematical complexity of viscoplastic fluids . . . . 51.2 VPL flows: Origin of Macro-size drop encapsulation . . . . . 91.3 Particles, bubbles and drops in viscoplastic fluids . . . . . . 111.4 Outline of thesis . . . . . . . . . . . . . . . . . . . . . . . . . 152 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . 162.1 Physical description . . . . . . . . . . . . . . . . . . . . . . . 162.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . 182.3 Dimensionless governing equations . . . . . . . . . . . . . . . 182.4 The plane Poiseuille solution UP (y) . . . . . . . . . . . . . . 212.5 Inertial effects and the droplet flow . . . . . . . . . . . . . . 212.6 Why large droplets yield the plug . . . . . . . . . . . . . . . 223 Asymptotic Solution . . . . . . . . . . . . . . . . . . . . . . . . 253.1 Rescaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2 Leading order solution in the yielded region . . . . . . . . . . 273.3 O(δ) solution in the yielded layer . . . . . . . . . . . . . . . 28ivTable of Contents3.4 Size of the yield surface perturbation from yy,0 and effects ofinertia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 304 Numerical Simulation . . . . . . . . . . . . . . . . . . . . . . . 334.1 Assumptions and simplifications. . . . . . . . . . . . . . . . . 334.2 Rheolef . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 354.3 Mesh adaptation . . . . . . . . . . . . . . . . . . . . . . . . . 354.4 Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.5 Effects of droplet spacing . . . . . . . . . . . . . . . . . . . . 395 Results and Discussion . . . . . . . . . . . . . . . . . . . . . . 415.1 Encapsulation and failure . . . . . . . . . . . . . . . . . . . . 415.2 Maximum size of encapsulated droplets . . . . . . . . . . . . 435.3 Encapsulation with fluids of different densities . . . . . . . . 476 Axisymmetric Geometry . . . . . . . . . . . . . . . . . . . . . 516.1 Poiseuille flow solution . . . . . . . . . . . . . . . . . . . . . 526.2 Slender drop . . . . . . . . . . . . . . . . . . . . . . . . . . . 526.3 Examples of encapsulation and failure . . . . . . . . . . . . . 556.4 Maximum size of encapsulated droplets . . . . . . . . . . . . 557 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 597.1 Summary of the results . . . . . . . . . . . . . . . . . . . . . 597.2 Limitations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 627.3 Future directions . . . . . . . . . . . . . . . . . . . . . . . . . 63Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64A Algebraic Details . . . . . . . . . . . . . . . . . . . . . . . . . . 71A.1 Channel Poiseuille flow of a Bingham fluid . . . . . . . . . . 71A.2 Asymptotic form of Bingham constitutive relations . . . . . 73A.2.1 channel geometry . . . . . . . . . . . . . . . . . . . . 73A.2.2 pipe geometry . . . . . . . . . . . . . . . . . . . . . . 74A.3 Perturbations are O(δ2) . . . . . . . . . . . . . . . . . . . . . 75A.3.1 yT − yy,0 ∼ O(δ2) . . . . . . . . . . . . . . . . . . . . 75A.3.2 U0 + δU1 = UP +O(δ2) . . . . . . . . . . . . . . . . . 76B Unsuccessful attempts . . . . . . . . . . . . . . . . . . . . . . . 78B.1 Method of islands of broken plug of [23] . . . . . . . . . . . 78B.2 Finding an admissible form of stress . . . . . . . . . . . . . . 80vList of Figures1.1 Examples and mathematical models of Viscoplastic fluids . . 41.2 Illustration of VPL flows . . . . . . . . . . . . . . . . . . . . . 101.3 “Frozen-Into” patterns in VPL flows . . . . . . . . . . . . . . 122.1 Schematic description of the problem . . . . . . . . . . . . . . 173.1 Example of asymptotic results, pressure and velocity andyield surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.2 Asymptotic results: yield surface . . . . . . . . . . . . . . . . 324.1 Schematic description of the problem . . . . . . . . . . . . . . 344.2 Mesh adaptation cycles . . . . . . . . . . . . . . . . . . . . . 374.3 Code validation . . . . . . . . . . . . . . . . . . . . . . . . . . 384.4 Effect of spacing . . . . . . . . . . . . . . . . . . . . . . . . . 394.5 Variation in the pressure drop with spacing between droplets 405.1 Speed distribution(u) as the droplet gets larger . . . . . . . . 425.2 Stress distribution as the droplet gets larger (slender drops) . 445.3 Stress distribution as the droplet gets larger (fat drops) . . . 455.4 Maximum size of encapsulated droplets . . . . . . . . . . . . 465.5 Effect of density difference on encapsulation . . . . . . . . . . 495.6 Stress distribution as the droplet becomes heavier . . . . . . . 506.1 Stress distribution in the axisymmetric geometry as the dropletgets longer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 566.2 Stress distribution in the axisymmetric geometry as the dropletheight is increased . . . . . . . . . . . . . . . . . . . . . . . . 576.3 Maximum size of encapsulated droplets axisymmetric geometry 587.1 Encapsulation of drops with exotic shapes. . . . . . . . . . . . 60B.1 Methods of island of broken plug . . . . . . . . . . . . . . . . 79viAcknowledgementsMy special thanks goes to my supervisor, professor Ian Frigaard for hisknowledge, his constant supports and patience, his caring and efforts. Itwas a privilege for me to work under his supervision. I cannot assume tohave any better scientific mentor. Thanks Ian!I would like to thank a fabulous couple, Drs. Sarah Hormozi and HamedGhasvari, who changed my life with their guidance and advice. Sarah as-sisted me in the analytical part of my project as a co-supervisor. She isalways a model of a hard-working, energetic and passionate researcher. Wehad many many insightful discussions with Ian and Sarah through the courseof this project.I also acknowledge kind helps of Ali Roustaei. Ali just simply knowseverything in the world of computation. Whenever I faced a problem, Iknew Ali has some hints. In addition, all my other officemates in ComplexFluids Lab, especially Ida without whom I do not have any idea how daysof graduate studies would look like, thank you all.Last but not the least, I express my gratitude to Dr. Green, the head ofmechanical engineering department at UBC.viiDedicationTo my ever-lasting friend, with whom I share childhood memories andgrown-up dreams, who has been the never-ending source of love and kind-ness, to my brother,To AminviiiChapter 1IntroductionEncapsulation is a process of entrapping one substance in another one. Themain function of encapsulation is to isolate an aggressive component fromthe environment and/or deliver it to a particular receptor [47]. This couldbe beneficial in transportation, coating, site-specific drug delivery, medicalimaging as well as food, cosmetic and pharmaceutical product manufactur-ing [29, 30, 12, 76].Various techniques of encapsulation are proposed in the literature; e.g.[43, 26, 12, 47, 39, 74, 73], but it is not our intention to review these here.A common feature of many current techniques is that they are limited todroplet sizes governed by capillary forces. This means that only small dropsize (order of millimetre or even smaller) can be encapsulated successfully.Instead, here we propose a method of encapsulation focused at Macro-scaledroplets, independent of capillary forces. The main advantage of this methodis that significantly larger droplets may be stably encapsulated, governed bythe length scale of the flow and yield stress of the encapsulating fluid. Unlikemany encapsulation processes, here we specially mean that the droplet isheld within the “unyielded plug region of a viscoplastic fluid”. These termsare defined later.This will open up use of encapsulation in many large scale industrial ap-plications. In particular, oil and gas, food, medical and cosmetic industriescan benefit from this. One interesting example is a process called Acidizingwhich will be explained here.Acidizing oil and gas reservoir: an application of Macro-size dropencapsulationProduction of an oil or gas well depends on the reservoir’s porosity andpermeability. Porosity quantifies how much void area exists in the rock for-mation and permeability measures the ability of the rocks to transmit oil orgas. The more porous and permeable the well, the more hydrocarbon canbe extracted. One of most common ways of enhancing oil production is tostimulate the reservoir by pumping acidic fluids down the well. This process11.1. Viscoplastic fluidscalled acidizing increases both porosity and permeability through openingnew cracks in the rock formation and by dissolving carbonate materials be-tween the sediment grains of the reservoir rocks. The process of acidizing hasbeen long practiced. The history of acid treatment of oil wells is discussed by[40]. One of the earliest studies on acidizing using hydrofluoric acid is theclassic study of [69]. With the development of computational techniques,attempts to model this process started around 1970’s; see the referencesreviewed in [46]. More recent studies involve multi-scale multi-dimensionalmodelling of the process [11, 24, 25, 63, 49].In spite of all these developments, there are still some crucial problems:• Hydrochloric acid (HCL) and hydrofluoric acid (HF) are the two mostcommon acids in the acidizing process. These materials are detrimen-tally corrosive for the casing (metallic wall of the well).• Acidic fluids should start reacting at a certain location down insidethe well, meaning that acid reactions should be retarded.As an example, [2] reported an unsuccessful acid job which induced signifi-cant damage to the oil well. Various mechanisms of inhibition and retarda-tion are being developed in order to address these issues; e.g. see [62, 58, 75].Potentially if we can encapsulate acid drops inside a friendly, neutral,non-corrosive fluid, we can pump the acid down the well. Where dropsof acid reach the desired location, we can break the encapsulation and letthem react with the rocks. Breaking may occur naturally, e.g. throughhigh shear as the drops enter the fracture, or could be controlled chemically.This automatically solves both issues mentioned earlier. We have eliminatedcontact of acidic fluids with the metallic wall of the casing; in addition, theoperator can decide when and where the reactions should start. Note thatwith the fluid volumes pumped in oil industry applications, use of micro-sizeencapsulation is out of question.Utilizing a specific category of non-Newtonian fluids, called viscoplasticfluids, we will suggest a new method of encapsulation and study the feasi-bility of this method in accommodating larger drop size. Here we introduceviscoplastic fluids:1.1 Viscoplastic fluidsIn Newtonian fluids shear stress τˆij and strain rate ˆ˙γij are linearly pro-portional. The constant of proportionality is called viscosity µˆ which is a21.1. Viscoplastic fluidsphysical property of the fluid.τˆij = µˆˆ˙γij (1.1)This model defines fluid as a material which does not resist shear stress.Once the fluid experiences a shear stress, it starts to flow. Assuming thetemperature is kept constant, the viscosity of a Newtonian fluid does notchange with time, shear rate or shear stress, and remains always constant.However there are many fluids in our daily life as well as in industrialand engineering applications for which viscosity does change. These fluidsare generally termed non-Newtonian fluids. The one important categorythat we are interested in is called “viscoplastic fluids”, or interchangeably,“yield stress fluids”. In contrast to Newtonian fluid, these materials do notdeform unless a given yield stress (τˆY ) is exceeded. In other words, thereis a minimum stress below which the fluid acts like a rigid solid, resistingthe stress. But once yield stress is exceeded, the fluid starts to deformand flow. Examples of these fluids are dairy products, liquid chocolate,drilling mud and lava; see Figure 1.1a. Many mathematical relations areproposed to model behaviour of viscoplastic fluids. Two common ones arethe “Bingham” and “Herschel-Bulkley” models defined as follows:1Bingham Modelτˆij =(µˆ+τˆYˆ˙γ)ˆ˙γij ⇐⇒ τˆ > τˆY (1.2a)ˆ˙γ = 0⇐⇒ τˆ ≤ τˆY (1.2b)Herschel-Bulkley Modelτˆij =(κˆˆ˙γn−1+τˆYˆ˙γ)ˆ˙γij ⇐⇒ τˆ > τˆY (1.3a)ˆ˙γ = 0⇐⇒ τˆ ≤ τˆY (1.3b)1All through the thesis, we adopt the convention of denoting all dimensional variableswith a “hat”, i.e. ·ˆ, and dimensionless variables without31.1. Viscoplastic fluidsFigure 1.1: a) Examples of visco-plastic fluids. Top row (left to right); paint,toothpaste, hagfish mucus, Middle row (left to right); dairy products, liquidchocolate, diaper cream, Bottom row (left to right); basaltic lavas, polymergels, drilling Fluids b) Relationship between shear rate (ˆ˙γ) and shear stress(τˆ) for Bingham fluids (solid line) and Herschel-Bulkley fluids (dashed line).Both pictures are taken from [32].41.1. Viscoplastic fluidsin whichˆ˙γij =∂uˆi∂xˆj+∂uˆj∂xˆi(1.4a)ˆ˙γ =√√√√123∑i,j=1|ˆ˙γ2ij | (1.4b)τˆ =√√√√123∑i,j=1|ˆ˙τ2ij | (1.4c)Here µˆ and τˆY are the plastic viscosity and yield stress. The power law indexis n and κ is referred to as the consistency, which corresponds to the plasticviscosity when n = 1. The Bingham model simply assumes a constant(tangent) viscosity when the fluid is yielded. However, in Herschel-Bulkleymodel the viscosity of the fluid after yielding is a function of shear rate.This models a shear thinning (n < 1) or shear thickening (n > 1) materialwith a yield stress; the case of shear thinning is shown in Figure 1.1b.1.1.1 Poiseuille flow of a viscoplastic fluidThe fully developed laminar flow of a fluid inside a long pipe or between twoinfinitely-wide plates is called Poiseuille flow. Poiseuille flow of a Newtonianincompressible fluid results in a parabolic profile for velocity which takesits maximum at the center line. Due to the symmetry, the shear stress iszero at the center line and is maximum on the wall. Therefore if we usea viscoplastic fluid instead, because of the existence of a yield stress, aninner layer arises at the center of channel or pipe within which stress isbelow yield stress. Therefore fluid inside the layer does not deform, and iscommonly called the “plug region”. Velocity within the plug region doesnot vary and the fluid resembles a moving solid object. We will analyticallysolve Poiseuille flow of a Bingham fluid inside a channel in §A. Mathematical complexity of viscoplastic fluidsFrom a mathematical point of view, modelling viscoplastic fluids is intrin-sically difficult. Here we explain what sort of complexities have to be dealtwith:51.1. Viscoplastic fluidsAnalytical difficultyThin film analysis and lubrication approximation have been widely utilizedin order to study problem involving motion of particles, drops or bubbles;e.g. see [13]. However, when dealing with viscoplastic fluids, it is known thata na¨ıve application of similar approximation methods can lead to solutionsthat exhibit velocity gradients in the plug region. This is inconsistent withthe Bingham model which defines plug regions as strain-rate-free domains.The inconsistency is called the “lubrication paradox” which has led to manyconfusing and incorrect statements. Confronted by this paradox, [42] con-cluded that no true plug region can exist in the flow of Bingham fluid incomplex geometries. They added that, in spite of approximately shear-freeplug regions, yielding must happen everywhere. This issue also questionedthe physical robustness of the Bingham fluid as mathematical model.Generally, the lubrication paradox appears as a result of a wrong scalinginside the plug region. In a long thin geometry, inside the yielded regionscales of pressure and shear stress are determined by balancing, to the lead-ing order, the shear force and pressure gradient. As we enter the plug,these scales are likely to remain valid. However, the extensional stresseswithin the plug become of the same order as the shear stresses. Therefore,in order to resolve the lubrication paradox, different scales must be intro-duced inside the plug. This led to development of more accurate methodswhich confirmed existence of true plug regions in complex geometries. Asexamples:• The earliest attempt to resolve the lubrication paradox was the workof [72] in which they studied the axial flow of a Bingham fluid inan eccentric annulus. They adopted a very systematic perturbationtechnique and cleverly defined a plug-like region where stress exceedsyield stress but only by a very small amount and called it “pseudo-plug”2. With a proper definition of scales in the pseudo-plug 3 theymanaged to show that true plugs exist at the widest part of the annulargap (as well as narrowest part).• [3] considered the case of gravity driven flow of a Bingham fluid over2The idea of plug-like region where stress just exceeds the yield stress was first usedby [53]3Inside the pseudo-plug velocity with the form of wpp = wpp0(θ) + δwpp1(δ, r) +O(δ2)is chosen where δ is the asymptotic parameter of the problem. This allows O(δ) strainrate inside the pseudo-plug. The pseudo-plug is an enclosed region which requires a thirdregion where the stress has to gradually relax from −B to B. This is where τ < B andthe true plug exists.61.1. Viscoplastic fluidsa slope. Incorporating a similar approach to [72], they argued thatno true plug could exist in such a geometry. However, there exists apseudo-plug region. In fact, the leading order velocity here is deter-mined by the film height and therefore, as the height of film changes,slow velocity gradients arise in flow direction. Thus, absence of thetrue plug is expected. The author concluded that existence of trueplug regions in varying geometries is more likely in wall-bounded ge-ometries.• [23] examined flow of a Bingham fluid between two wavy plates. Theirgoal was to study channel Poiseuille flow where small perturbationsare introduced on the wall of the channel; in particular, how the plugregion changes in response to the perturbations. First of all, theyshowed that the conventional lubrication scaling ends up with the lu-brication paradox (variation of velocity in principal direction withinthe plug). To resolve this inconsistency, two different scales are pro-posed inside the pseudo-plug depending on the relative configurationof the pseudo-plug and the true plug and finally the existence of trueplug region is confirmed provided that the amplitude of perturbationsis sufficiently small. With the aid of the proper scaling inside theplug, they also estimated a sufficient condition for the plug to yield.More recently, [61] revisited this problem and verified the asymptoticsolution provided by [23].Computational difficultyViscoplastic models are all singular when strain rate vanishes, which im-plies an infinite effective viscosity in unyielded flow regions. In addition, theposition of plug interface is not known in priori. To deal with the singular-ity problem two different classes of methods are proposed: Regularisationmethods and augmented Lagrangian algorithms. In general, viscosity regu-larization methods suffer from an uncertainty, especially in determining theplug interface, while augmented Lagrangian methods are more complicatedin implementation and more costly computationally [15, 4].• Class One, Regularisation Methods: In regularisation methods, weeliminate the singularity of stress by replacing the plug region with ahighly viscous fluid (as opposed to the original model which considersthe plug region a solid). With this assumption we get rid of the singularpoint (with infinite derivative) near origin, however, it still remainsnon-differentiable. The efficiency of different regularisation methods71.1. Viscoplastic fluidsare compared in [22]. As an example, the “simple” regularisationmodel is written as:τˆij =(µˆ+τˆYˆ˙γ1 + ˆ˙γ0)ˆ˙γij (1.5)where ˆ˙γ0 is a representative strain rate of the flow and in which  1. Thus, the regularisation is negligible if the strain rate is oforder ˆ˙γ0. As the strain rate approaches zero, the effective viscosityapproaches τˆY /(ˆ˙γ0), remaining large but finite. Compared to Aug-mented Lagrangian method (introduced below), regularisation meth-ods are computationally cheaper and easier to implement [15]. How-ever, as  −→ 0, the method becomes costly. Unfortunately in spiteof many non trivial problems computed with this method, there isno theory which guarantees regularisation methods converge to thecorrect stress field.In addition, regularisation methods may fail to give the correct posi-tion of the yield surfaces, and potential errors are largest in flows withlong-thin geometries (as is our case); see [22]. This is probably themain drawback of regularisation methods [61]. For flows such as thathere, where we investigate whether or not the droplet is fully encapsu-lated in unyielded fluid, regularization methods are unable to answerwith certainty.• Class Two, Augmented Lagrangian (AL) Algorithms: The alterna-tive to regularization methods is algorithms such as the augmentedLagrangian method of [20, 27]. Through the augmented Lagrangianalgorithms, the exact form of constitutive relation is implemented inthe discretized equations which gives the true unyielded regions withexactly zero strain rate. As we want to study the perturbations ofplug interface, and break of the plug, it is crucial to distinguish thetrue unyielded regions of the flow from those where the strain rate ismerely small.As an example, let us consider the augmented Lagranian method whichis developed for Bingham fluids. This method is based on the vari-ational form of Stokes equation for Bingham plastic, derived e.g. by[19], which says that the velocity field of a Bingham fluid minimisesthe following variational form:Jˆ(vˆ) =µˆ2∫Ωˆ||ˆ˙γ(vˆ)||2 + τˆY∫Ωˆ||ˆ˙γ(vˆ)|| −∫Ωˆgˆvˆ (1.6)81.2. VPL flows: Origin of Macro-size drop encapsulationin which ˆ˙γ(vˆ) = ∇vˆ and ||ˆ˙γ(vˆ)|| is the norm of ˆ˙γ(vˆ) and gˆ repre-sents a body force. Ωˆ is also the domain of the problem. The non-differentiable term of∫Ωˆ ||ˆ˙γ(vˆ)|| hampers the use of the standard opti-mization techniques. Instead, two Lagrange multipliers (γˆ, Tˆ) can beintroduced to relax the problem into a saddle-point problem [20]:Jˆa(vˆ, γˆ, Tˆ) = Jˆ(vˆ) +∫Ωˆ(ˆ˙γ(vˆ)− γˆ) : Tˆ +aˆ2∫Ωˆ||ˆ˙γ − γˆ|| (1.7)where aˆ is the augmented Lagrangian parameter. To find the saddlepoint, a Uzawa-type algorithm similar to the one in [67] is useful. Ithas been shown that T and γ converge to “admissible” stress andstrain tensors. A modern review of such methods is given by [28].1.2 VPL flows: Origin of Macro-size dropencapsulationThe idea of Macro-size drop encapsulation stems from observations in aseries of multi-layer flow studies called VPL flows. Viscoplastically Lubri-cated (VPL) flows are a specific arrangement of multi-layer flows where ayield stress fluid is used to stabilize the interface of the two layers. It isknown that interfacial instabilities inherently occur in multi-layer flows asthe physical properties of the fluids change across the interface. If we canposition a visco-plastic fluid in the outer layers and tune the flow rates tokeep the interfacial shear stress less than the yield stress (τY ), the fluid-fluidinterface of the flow is converted to a fluid-solid interface, which does notallow instabilities to develop. In other words, a visco-plastic fluid adjacentto the wall lubricates a centrally positioned viscous fluid (see Figure 1.2).This idea was originally introduced in [21] in which the author consid-ered a two-layer flow of two Bingham fluids where the inner fluid has asmaller or zero yield stress. Through stability analysis, he demonstratedthat this configuration of multi-layer flow does not suffer from the classicallinear interfacial instabilities that are found for analogous multi-layer flowsof Newtonian and many other non-Newtonian fluids. Instead, such flows firstbecome linearly unstable to shear instabilities that are analogous to thosefound for single fluid flows. This is then continued by [51, 52, 36, 34] wherethey studied the nonlinear stability of a multi-layer shear flow, for New-tonian and non-Newtonian core fluid, and demonstrated that with carefulchoice of rheology, multi-layer flows can be achieved which are linearly andnon-linearly stable even at significantly high Reynolds numbers and large91.2. VPL flows: Origin of Macro-size drop encapsulationFigure 1.2: Illustration of VPL flows. a) Visco-plastically lubricated pipeflow; b) Cross-sectional view c) Schematics of the type of basic velocityprofiles considered, with reduced rate of strain in the Newtonian fluid andan unyielded plug zone adjacent to the interface in the lubricating fluid; [32]amplitudes of perturbation. The transient characteristics (entry, start upand development length after which the flow is parallel) of these flows arealso studied in [37]. More interestingly, [36] computationally confirms thatstable VPL flows can be achieved with more layers (e.g. seven layers). Inaddition, they have experimented numerically with different injection con-figurations in order to advect encapsulated “letters” and “write” in the flow.This is where a vague idea of encapsulation is introduced, but not treatedvery systematically.[38] studied the feasibility of VPL flows experimentally. The schematicsof their experimental setup is shown in Figure 1.3a. Xanthan4 as the innerlayer and Carbopol as the outer layer5 flow in path 1 and path 2 respectively.Later [35] replaced the core fluid with Polyethylene oxide (PEO) which ex-hibits elastic behaviour to show that the VPL technique is also applicablefor multi-layer flows with elasticity properties. The idea of VPL flows hasbeen extended in different ways. For example, recently [18] and [33] haveintroduced periodic perturbations to the flow rates of initially stable VPLflows, resulting in stable patterned interfaces. The setup shown in Figure4Xanthan has shear-thinning rheology.5Carbopol is the most common viscoplastic fluid in experimental practices.101.3. Particles, bubbles and drops in viscoplastic fluids1.3a has changed such that a temporal profile of flow rate can be givento the inner or outer layer. Any sudden change in the flow rate producessome instabilities in the flow which could form frozen-into shapes within theinterface of VPL flow. For example, Figure 1.3 demonstrates one typicalexperiment of [18] where a “diamond necklace” pattern is preserved.We will extend this idea by considering the encapsulation of dropletswithin a visco-plastic fluid.1.3 Particles, bubbles and drops in viscoplasticfluidsAlthough it is not directly relevant to our study, it is still useful to reviewthe previous works of motion of particles, bubbles and drops in viscoplasticfluids. Motion of particles, bubbles and drops in viscoplastic fluids is becom-ing an ubiquitous process in many engineering and industrial applicationssuch as oil and petroleum, pharmaceutical, cosmetic and medical industries.For instance:1. Drilling mud used in the construction of oil and gas wells is considereda viscoplastic fluid. If the hydrostatic pressure of the mud is unex-pectedly lower than the pore pressure of gas in the surrounding rockformation, gas enters and rises upward in the well. This phenomenonis called gas kick and can eventually cause a blowout at the surfacewith enormous environmental and economical costs.2. In the processing and packaging stages of various products bubblesmay either be unwanted or desirable, e.g., hair gel is sold by vol-ume (bubbles included). In contrast, it is not always desirable tohave air trapped in chocolate as aerated chocolate has a different tastethan solid chocolate. Also bubbles may affect mechanical propertiesof structural materials like cement. [16, 15].3. In materials such as cement or toothpaste which contain suspendedparticles, sedimentation breaks the homogeneity of the fluid and there-fore alters the physical properties [54].In the absence of fluid motions there are many studies of settling in sta-tionary yield stress fluid whereby the yield stress holds a particle stationary111.3. Particles, bubbles and drops in viscoplastic fluidsFigure 1.3: a) Schematic view of [18] experimental setup. b) Formation ofa stable “diamond necklace” pattern in his experiments. Inner flow: PEOis initially flowing at a rate of Qinner = 15.4ml/s, suddenly the flow rateis reduced to Q′inner = 2.6ml/s, held at this flow rate for T = 1s and thenreturned back to its initial value. These pattern repeats every 10 seconds.Outer flow: Carbopol flowing at the rate of Qouter = 8.6ml/s.121.3. Particles, bubbles and drops in viscoplastic fluidsuntil a certain critical buoyancy stress is exceeded. Thus, particles below acertain size are held in suspension. In the case of a spherical particle, [6]simulated the flow around the particle using a regularization method andinterestingly, found out that the spherical particle falls only ifYg =τˆYRˆgˆ∆ρˆ≤ 0.095. (1.8)where Rˆ is the radius of sphere and ∆ρˆ and gˆ are density difference and grav-itational acceleration respectively. Their results are later confirmed by manyother studies showing only minor discrepancies; e.g. [8, 5, 44]. More recently,[60] analyzed Stokes flow around a particle in Bingham fluid with emphasison the resistance and motility problem. Later,[59] performed careful ex-periment on the sedimentation of spherical particles. The main conclusionthey draw was inaccuracy of common viscoplastic models in characterizingthe rheological behaviour of viscoplastic material (in particular, Carbopol)which may lead to discrepancies between experimental and numerical data.The motility bound given in (1.8) does not hold for a flowing yield stressfluid [54]. This arises another class of studies which consider particle settlinginduced by shear. [50] claimed that in pipe flows no settling occurs in asheared yield stress fluid if Yg ≥ 6. However, based on their observations,[54] raised doubts as to the validity of this stability bound.As for cylindrical particles, in a theoretical paper, [1] approximated theBingham flow around a circular cylinder and indicated the main character-istics of velocity and stress fields. This problem has been later revisitedby [14] and [71] who utilized regurisation methods to simulate the problem.The later carefully studied the effect of different numerical parameters onthe convergence characteristics. In addition, they gave estimates for thedrag coefficient of the particle as a function of Bingham number. With theaid of augmented Lagrangian method [64] computed this problem when thecylindrical particle is confined by two parallel plates. They investigated theflow characteristics, in particular shapes of yielded and unyielded regions,when the particle is asymptotically close to the wall. The results of thispaper can be useful for the purpose of benchmarking and code validation.Studies on the rise of bubbles in a viscoplastic medium are not unfor-tunately as clear. The mathematical analysis of this category of flow wasinitiated by [7] and then continued in [70]. They both employed pertur-bation analyses and provided solutions for a bubble moving in a Binghamfluid for small Bingham numbers. Utilizing variational analysis, [16] deriveda sufficient condition (not necessary) for stopping of bubbles in Herschel131.3. Particles, bubbles and drops in viscoplastic fluidsBulkley fluids. Later, they [17] experimentally confirmed their analyticalbounds and proposed less conservative empirical limits.[68] simulated buoyancy driven flow induced inside a stagnant pool ofBingham fluid by a 2D bubble or drop. Using a regularization method,they computed the shape of plug region and argued the presence of rotatingunyielded zones at the equator of the bubble. Accuracy of regularizationmethods in simulating bubble rise in a viscoplastic fluid has been studiedby [15] who compared the regularization method introduced by [56] withthe results of the augmented Lagrangian method in Herschel-Bulkey model.They showed that up to intermediate Bingham numbers the regularizationmodel of [56] predicts the unyielded region accurately. However, for largerBingham number extra caution is required in choosing the regularisationparameter.In so far as viscous droplets are concerned, there are fewer studies inthe literature. For instance, [57] simulated motion and deformation of asingle drop in a Bingham fluid under gravity. They reported approachinga quasi-steady state after a relatively large transient period. The transientperiod and velocity magnitude depend strongly on Bingham number. Theyalso observed stationary flows for sufficiently large Bingham numbers. [31]experimentally studied the fall of a Newtonian drop inside an otherwisestagnant Bingham fluid. Using a combined PTV and PIV technique, theydetermined the approximate position of plug interface as well as the walleffect on sedimentation velocity.Interaction of multiple drops or bubbles in a viscoplastic media has alsorecently received attention [45, 57, 68, 41, 31]. In particular, [45] reportedthat the travelling particles in a viscoplastic media exhibit negligible inter-action unless they are in close proximity (less that four sphere radii). [57]examined the interaction of two drops while they are falling under gravity ina Bingham fluid. For the case of two similar drops, they showed that withinthe proximity range defined in [45], the drops tend to approach each otherand coalesce. More interestingly, [68] found in their simulations that uponcertain arrangement, a group of bubbles can rise under conditions where asingle drop is unable to move. In the context of our study, the proximitydistance of [45] is relevant in that we also will determine a droplet spacingabove which the droplets do not influence one another. However, note thatour droplets do not yield the fluid encapsulating fluid, so the connection withthe above studies is via locality of the stress field rather than the velocityfield.Finally, moving to large numbers of particles, bubbles or droplets weenter the realm of yield stress suspensions, foams and emulsions. Recent141.4. Outline of thesiswork has addressed the question of estimating the bulk yield stress of ayield stress suspension, according to the fraction of the dispersed phase.[10] developed a theoretical model to estimate the overall yield stress of asuspension of noncolloidal and non-Brownian particles immersed in a yieldstress fluid. They showed that such a suspension can be modelled as “aHerschel-Bulkley fluid with an exponent equal to that of the suspendingfluid.” This agrees reasonably well with the experimental results of [55, 48].Of course, in dealing with suspensions the dispersed phase length-scale ismuch smaller than that of the bulk flow, unlike here where the droplets areof the flow dimension and we focus on a continuous stream.1.4 Outline of thesisAn outline of the thesis is as follows. In chapter 2 we describe the problemof study and derive the governing equations in dimensionless form. Forsimplicity we start with the channel geometry. Chapter 3 considers slenderdroplets, showing that the yield surface is barely perturbed by the presenceof the droplets. In chapter 4 we describe the numerical approach that wetook to compute the flow around an iso-dense droplet. Later, in chapter 5 wepresent the results of computational analysis. In particular, we determinethe mechanisms of yielding as droplets grow large and find the maximaldroplet size. We also explore the effect of density difference on encapsulationprocess. Next, in chapter 6 we generalize the foregoing results for planechannel encapsulation to the more practical pipe geometry. The thesis closeswith a summary of contributions and future directions in chapter 7.15Chapter 2Problem StatementHere we describe the problem and specify the effective parameters. Wesimplify the problem by making some assumptions which help us solve itanalytically and numerically.2.1 Physical descriptionWe consider an infinite number of drops of evenly spaced flowing within ayield stress fluid inside an infinitely vertically-oriented long channel. Thewidth of channel is 2Rˆ and length and width of drops are 2Hˆ and 2Lˆ re-spectively. Aspect ratio of drops is defined as δ = HˆLˆ; see Figure 2.1 for aschematic view of the problem. The drops are to be considered Newtonianand the carrier fluid is modelled as a Bingham fluid (simplest model for aviscoplastic fluid). Both the bulk (carrier fluid) and disperse (drop) phaseare incompressible and the areal flow is 2RˆUˆ0, i.e. the mean velocity alongthe channel is Uˆ0. The droplets are spaced a distance lˆ apart, each centredat xˆ = klˆ + UˆP tˆ : k ∈ Z, and have boundaries (interfaces) given by:yˆ = ±hˆ(xˆ− klˆ − UˆP tˆ), (2.1)where hˆ(0) = Hˆ.As explained before in §1.1, in the absence of drops, velocity of the carrierfluid adopts a uniform symmetric Poiseuille profile, which in the case ofBingham fluid, contains an inner plug region in the middle of channel ofwidth 2yˆy0 and speed of UˆP where shear rate is zero and the fluid is justlike a moving solid. Insertion of drop inside the plug region will changethe stress distribution within the plug region, and possibly may cause it toyield (total stress exceeds the yield stress τˆY ). We want to see when thisplug region may encapsulate a fluid drop without the plug yielding. In thisthesis, when we refer to encapsulation, we specifically mean that the dropletis held in an unyielded plug region.162.1. Physical descriptionxˆRˆ2yˆLˆ2( )yU ˆˆ0,ˆ2 yyRlˆgˆHˆ2PlugYielded layersFigure 2.1: Schematic description of the problem172.2. Governing equations2.2 Governing equationsInside the drop continuity and momentum equations write:0 =∂uˆ∂xˆ+∂vˆ∂yˆ(2.2a)ρˆd(∂uˆ∂tˆ+ uˆ∂uˆ∂xˆ+ vˆ∂uˆ∂yˆ)= −∂pˆ∂xˆ+ µˆd∂2uˆ∂xˆ2+ µˆd∂2uˆ∂yˆ2+ ρˆdgˆ (2.2b)ρˆd(∂vˆ∂tˆ+ uˆ∂vˆ∂xˆ+ vˆ∂vˆ∂yˆ)= −∂pˆ∂yˆ+ µˆd∂2vˆ∂xˆ2+ µˆd∂2vˆ∂yˆ2(2.2c)where ρˆd and µˆd are the drop’s density and viscosity. Similarly in the carrierfluid0 =∂uˆ∂xˆ+∂vˆ∂yˆ(2.3a)ρˆ(∂uˆ∂tˆ+ uˆ∂uˆ∂xˆ+ vˆ∂uˆ∂yˆ)= −∂pˆ∂xˆ+∂τˆxx∂xˆ+∂τˆxy∂yˆ+ ρˆgˆ (2.3b)ρˆ(∂vˆ∂tˆ+ uˆ∂vˆ∂xˆ+ vˆ∂vˆ∂yˆ)= −∂pˆ∂yˆ+∂τˆxy∂xˆ+∂τˆyy∂yˆ(2.3c)where ρˆ is the carrier fluid’s density. In order to close (2.3) we need toincorporate the constitutive relations for a viscoplastic fluid. Taking carrierfluid as a Bingham fluid, the constitutive relation is (1.2) in which ˆ˙γ =√ˆ˙γ2xy + ˆ˙γ2xx, τˆ =√τˆ2xy + τˆ2xx, andˆ˙γxy = ˆ˙γyx =∂uˆ∂yˆ+∂vˆ∂xˆ(2.4a)ˆ˙γxx = −ˆ˙γyy = 2∂uˆ∂xˆ= −2∂vˆ∂yˆ(2.4b)To simplify these sets of equations and find the effective dimensionless pa-rameters, we introduce a scale for each parameter and nondimensionalizethe governing equations. As expressions of boundary conditions in the di-mensionless form look simpler, we explain them in the next section.2.3 Dimensionless governing equationsWe consider flows, periodic in xˆ and symmetric with respect to the chan-nel centerline yˆ = 0. We scale all lengths with Rˆ, velocities with Uˆ0 and182.3. Dimensionless governing equationsdeviatoric stresses with µˆUˆ0/Rˆ. The pressure is scaled as follows:pˆ = pˆ0 + ρˆgˆxˆ+µˆUˆ0Rˆp (2.5)i.e. p is the dimensionless modified pressure, subtracting off the static pres-sure ρˆgˆxˆ and a possibly time dependent pressure pˆ0(tˆ).It is assumed that the distance between droplets lˆ is large enough for thecarrier fluid to assume its undisturbed uniform Poiseuille profile betweendroplets, in which the plug moves with speed UˆP . We explore effects ofvarying lˆ later in §4.5. To fully exploit the symmetry of the problem, wetranslate to a moving frame as follows:(x, y) =(xˆ− UˆP tˆRˆ,yˆRˆ)(2.6a)t =UˆP tˆRˆ(2.6b)(u, v) =(uˆ− UˆPUˆ0,vˆUˆ0)(2.6c)In the moving frame, we consider the flow to be steady and due toperiodicity consider only the domain: (x, y) ∈ [−l/2, l/2] × [−1, 1], wherel = lˆ/Rˆ. The droplet-encapsulating fluid interface is given by: y = ±h(x)with h(0) = H = Hˆ/Rˆ and where h(x) = 0 for |x| ≥ L = H/δ = Lˆ/Rˆ.Adopting these scales, the dimensionless form of (2.2) writes:0 =∂u∂x+∂v∂y(2.7a)φRe(u∂u∂x+ v∂u∂y)= −∂p∂x+m∂2u∂x2+m∂2u∂y2+ χ (2.7b)φRe(u∂v∂x+ v∂v∂y)= −∂p∂y+m∂2v∂x2+m∂2v∂y2(2.7c)here φ = ρˆd/ρˆ is the density ratio, m = µˆd/µˆ is a viscosity ratio and χ is adimensionless group:χ =(ρˆd − ρˆ)gˆRˆ2µˆUˆ0. (2.8)192.3. Dimensionless governing equationsrepresenting the ratio of buoyant stresses to viscous stresses. The scaledequations of motion in the carrier fluid are:0 =∂u∂x+∂v∂y(2.9a)Re(u∂u∂x+ v∂u∂y)= −∂p∂x+∂τxx∂x+∂τxy∂y(2.9b)Re(u∂v∂x+ v∂v∂y)= −∂p∂y+∂τxy∂x+∂τyy∂y(2.9c)where the constitutive relations are:τij =(1 +Bγ˙)γ˙ij ⇐⇒ τ > B (2.10a)γ˙ = 0⇐⇒ τ ≤ B (2.10b)in which γ˙ =√γ˙2xy + γ˙2xx, τ =√τ2xy + τ2xx, andγ˙xy = γ˙yx =∂u∂y+∂v∂x(2.11a)γ˙xx = −γ˙yy = 2∂u∂x= −2∂v∂y(2.11b)in which all shear rates are scaled by Uˆ0/Rˆ. Two dimensionless groupsintroduced above are defined as: Re = ρˆUˆ0Rˆ/µˆ is the Reynolds number andB = τˆY Rˆ/µˆUˆ0 is the Bingham number.At the walls of the channel no-slip conditions are to be satisfied and atx = ±l/2 the flow is fully developed and adopts the Poiseuille profile:u(x,±1) = −up,0, v(x,±1) = 0, (2.12a)u(±l/2, y) = UP (y)− up,0, v(±l/2, y) = 0, (2.12b)Here UP (y) denotes the fully developed plane Poiseuille profile, which hasdimensionless plug speed up,0; see §2.4. At the interface between dropletand carrier fluid the velocity and the traction are continuous:JτntK = 0, at y = ±h(x), (2.13a)J−p+ τnnK = 0, at y = ±h(x). (2.13b)Here the subscripts n and t on the deviatoric stress components refer todirections resolved normal and tangential to the interface. For simplicity weneglect capillary effects in the formulation as the aim is to study encapsu-lation via yield stress effects on a macro-scale, i.e. implicitly meaning thatthe surface tension is insignificant in comparison to the yield stress.202.4. The plane Poiseuille solution UP (y)2.4 The plane Poiseuille solution UP (y)In §A.1 we straightforwardly show that channel Poiseuille flow solution of aBingham fluid can be resolved as:UP (y) =up,0 if |y| ≤ yy,0,up,0[1−(|y| − yy,0)2(1− yy,0)2]if |y| > yy,0,(2.14)in which up,0 is the plug velocity and is given byup,0 =B2yy,0(1− yy,0)2, (2.15)and yy,0, the position of the yield surface, is found by this cubic equationknown as Buckingham Equation:y3y,0 − (3 +6B)yy,0 + 2 = 0, (2.16)provided that the unity flow rate is satisfied. It has been shown in §A.1 that2.16 has a single root yy,0(B) ∈ (0, 1), and exhibits the following asymptoticbehaviour:yy,0 ∼B3−B26as B → 0; yy,0 ∼ 1−√2B1/2+O(B−1) as B →∞(2.17)2.5 Inertial effects and the droplet flowAs discussed previously, the aim of our study is to understand situationsunder which there may exist a droplet fully encapsulated in the unyieldedplug of the encapsulating fluid. What we consider below neglects inertialeffects. We now examine this assumption.Firstly with respect to the droplet flow, for any droplet that is fullyencapsulated within the plug the velocity at the interface satifies (u, v) = 0,as we have subtracted off the plug velocity. If we consider these conditionsas the boundary conditions for the droplet domain the droplet flow problemeffectively decouples from that of the carrier fluid. Using these conditions,the droplet flow has the (unique) Stokes flow solution (u, v) = 0, everywherewithin the droplet, withp = χx+ constant. (2.18)212.6. Why large droplets yield the plugThe trivial solution also solves the steady inertial problem. Assuming thatthe plug remains unyielded we may even consider time dependent solutionsfor the droplet. As there is no driving force for the flow within the droplet,an energy stability analysis would show that any transients in the dropletdecay to zero like ∼ exp−(φRe/m)t, which is the usual viscous decay timescaleof the droplet. There is no bound needed on the size of the initial conditionfrom the perspective of the decay bound. However, transients within thedroplet would induce stresses within the plug and hence a priori assump-tion of an unyielded plug could fail for sufficiently large initial conditions.Nevertheless, it is apparent that the trivial solution is likely to be stable fora reasonable range of Re > 0 and we therefore consider this to represent thedroplet solution.Since the fluid in the droplet is Newtonian, the trivial solution (u, v) = 0,implies vanishing of the shear stresses within the droplet. Equation (2.13)therefore simplifies to the following stress conditions to be satisfied by theBingham fluid at the interface:τnt = 0, at y = ±h(x), (2.19a)−p+ τnn = −χx+ constant, at y = ±h(x). (2.19b)Assuming symmetry and locating the droplet at the origin, we may considerthe constant on the right-hand side of (2.19b) to be zero.For the carrier fluid in the absence of the droplet, the Poiseuille flowsolution is valid up to large Re, when the flow destabilizes due to turbulenttransition. Departures from this solution in the droplet flow solution arelikely to be associated with perturbation of the streamlines in the vicinity ofthe droplet. Conditions (2.19) are satisfied at the interface, which results in aperturbation of the stress field within the plug. This in turn may perturb theyield surface from its uniform value yy,0 and the perturbed yield surface maythen induce a velocity perturbation from the Poiseuille flow solution. If wesuppose that the aspect ratio of the streamline perturbation is characterisedby a parameter ε then the inertial terms in the x-momentum equation havesize εRe and those in the y-momentum equation have size ε3Re. We mayconsider the inertial terms to be small if εRe  1, which we show later tobe the case.2.6 Why large droplets yield the plugIt might at first appear that the plug may sustain any droplet of size h(x) <yy,0, since the rigid motion of the uniform plug around the droplet allows222.6. Why large droplets yield the plugthe droplet to translate at uniform speed along the channel. However, thisreasoning neglects the effects of the droplet on the stress field. We consider2 simple examples that suggest how the stresses may act to break the plugregion for sufficiently large droplets.First let us consider an iso-dense droplet (χ = 0; see later in §5.3 for χ 6=0). We have seen that the undisturbed flow (no droplets) is the Poiseuilleflow of §2.4 above, in which the frictional pressure gradient satisfies:∂p∂x= −Byy,0, (2.20)Now let us introduce a droplet within the plug region, at x = 0 in thetranslating frame. If the droplet translates uniformly, the shear stresses inthe droplet are zero and the pressure in the droplet is simply p = constant,(chosen to be zero). In contrast, within the yielded layer of the carrier fluidwe have p = −Bx/yy,0, varying by ∼ 2BL/yy,0 along the length of thedroplet.The stresses within the unyielded plug are indeterminate in a yield stressfluid. Stress continuity implies that the tangential stress must vanish atthe interface, whereas the total normal stress must balance the constantdroplet pressure. We see that the pressure imbalance between the dropletand the yielded layer of the encapsulating fluid is of size ∼ BL/yy,0, whichmust somehow be absorbed by the deviatoric stresses within the plug region.Consequently, a simple order of magnitude estimate for the length of dropletthat can be sustained by a visco-plastic fluid is simply:yy,0 & L. (2.21)Putting this into dimensional terms we have:τˆYGˆRˆ&LˆRˆ, ⇒τˆYGˆ& Lˆ, (2.22)where Gˆ is the frictional pressure gradient of the uniform Poiseuille flow.Since yy,0 < 1 we see that (2.21) is fundamentally a restriction on dropletarea (equivalently volume in three dimensions), i.e. we expect that LH <yy,0.A second example considers the tangential stress distribution. Considera droplet with slowly varying height (∣∣∂h∂x∣∣  1), encapsulated inside theplug region. Inside the plug the stress distribution is undetermined, but hasto satisfy the momentum equations in (2.9). Assuming the yield surface is232.6. Why large droplets yield the plugapproximately yy,0, we might linearly approximate:τxy ≈ −By − hyy,0 − h, (2.23)i.e. satisfying the stress conditions at the droplet and yield surface. Ignoringthe x-variation of h, from (2.9b) we find that p−τxx = −Bx/(yy,0−h) withinthe plug. Outside the plug we have p = −Bx/yy,0, and imposing continuityof normal stress across the yield surface:−Bx/yy,0 = p− τyy = p+ τxx, ⇒ τxx =Bx2hyy,0(yy,0 − h). (2.24)Now observe that for fixed x, as yy,0 − h → 0, |τxx| → ∞, i.e. suggestingthat droplet widths may not approach arbitrarily close to the plug width.Finally if we assume the distributions (2.23) and (2.24), on combiningwith the yield criterion (τ = B) we predict that the plug yields along y =±yy(x):yy(x)− hyy,0 − h≈√1−x24( 1yy,0 − h−1yy,0)2(2.25)Thus, also x is limited in such a stress distribution, i.e. |L| < 2yy,0(yy,0 −h)/h, in order for the plug to remain unyielded around the droplet, or LH <2yy,0(yy,0−h). These two simple examples represent extremes of behaviour,indicating simply that we may expect limits on the size of encapsulateddroplets. In practice the stress distributions within the unyielded plug canbe any that satisfy the Stokes equations and boundary conditions. The ideaof constructing admissible form of stress distributions inside the plug regionis elaborated in more details in §B.24Chapter 3Asymptotic SolutionA common class of flows for which it is possible to construct analyticalsolutions is that in which the streamlines are nearly one-dimensional. In theproblem considered, non-uniformity only arises from the droplet shape andhence we look at droplets of small aspect ratio (H/L = δ  1). However,the analysis above leading to (2.21) suggests that (isodense) droplets cannot be encapsulated if L is larger than yy,0. Combined with the requirementof δ  1 implies that H ∼ δ  1, i.e. we consider small-slender droplets.Assuming h(x) ∼ O(δ), we develop an asymptotic approximation to theflow in the case that δ  1. Evidently we expect that the axial velocitywill have similar size to that of the droplet free flow and the appropriatescale for y remains the channel half-width. On the other hand, since theshear stress is zero at y = h(x), this suggests an O(δ) perturbation of thestress field, which may induce velocities v ∼ δ. In order to balance the massconservation equation, a rescaling of x is needed. Therefore, we consider thefollowing rescaling:3.1 Rescalingxδ = X, y = Y, u = U, v = δV, pδ = P. (3.1)The deviatoric stresses are rescaled according to the velocity scales andimplies:γ˙xy = γ˙XY , γ˙xx = δγ˙XX , τXY = τxy, τXX = δτxx (3.2)Note that this type of scaling leads to approximation based on the yieldedflow region and driven by the flow geometry and as it has been discussedearlier in §1.1.2 can lead to the lubrication paradox. Ignoring the inertial253.1. Rescalingterms, the governing equations after rescaling (3.1) write:0 =∂U∂X+∂V∂Y(3.3a)0 = −∂P∂X+ δ2∂τXX∂X+∂τXY∂Y(3.3b)0 = −∂P∂Y+ δ2∂τXY∂X+ δ2∂τY Y∂Y(3.3c)The constitutive relations are similar to (2.10) except with τ =√τ2XY + δ2τ2XXand γ˙ =√γ˙2XY + δ2γ˙2XX :γ˙XX = 2∂U∂X, γ˙XY =∂U∂Y+ δ2∂V∂X.Adopting a regular perturbation expansions in δ of form:(P,U, V ) = (P0, U0, V0) + δ(P1, U1, V1) + δ2(P2, U2, V2) + ...it has been shown in Appendix A.2.1 that the asymptotic form of Binghamconstitutive relations is resolved by:τXY =[∂U0∂Y+B sgn(∂U0∂Y)]+[δ∂U1∂Y]+O(δ2) (3.4a)τXX = 2(1 +B|∂U0∂Y |)∂U0∂X+O(δ) (3.4b)provided that the Bingham fluid is yielded. Assuming the above scal-ing, our periodic domain becomes X ∈ [−lδ, lδ], Y ∈ [−1, 1], althoughwe may consider only Y ≥ 0 through symmetry. The droplet occupiesX ∈ [−Lδ, Lδ] = [−H,H] and Y ∈ [−h(X), h(X)). No slip boundary con-ditions are satisfied at the walls. In parts of the channel where there is adroplet, the two conditions (2.19a) & (2.19b) are satisfied at Y = ±h(X).In terms of the rescaled variables these become:τXY − δ2 2τXXhX1 + δ2h2X= 0, (3.5a)P + δ2(1− δ2h2X)τXX + 2τXY hX1 + δ2h2X= χX, (3.5b)where hX is the derivative of h(X) with respect to X. Ignoring the tip ofthe droplet X → ±H, we may assume that |hX | ∼ O(1). Assuming that263.2. Leading order solution in the yielded regionthe droplet remains encapsulated in the plug region of the flow the stressesremain indeterminate within the plug, away from the interface. In order toapply these conditions to deriving a perturbation solution we need at leastto understand the order of magnitude of the stress components. The scalesused are derived from the yielded region, where they are determined by thevelocity scales and a leading order shear flow balance for the pressure. Aswe enter the unyielded plug, the scaling for P and τXY is likely to remainvalid. However as previously discussed in § 1.1.2, in long-thin geometriesthe plug yields the extensional stresses within the plug become of the sameorder as the shear stresses, i.e. this is the yielding mechanism. Therefore,although the leading order terms in (3.5a) & (3.5b) are clear:τXY,0 = 0 at Y = ±h(X), (3.6a)P0 = χX at Y = ±h(X), (3.6b)at first order there may be some adjustment close to breaking of the plug:τXY,1 − 2δτXX,0 = 0 at Y = ±h(X), (3.7a)P1 + δτXX,0 = 0 at Y = ±h(X). (3.7b)3.2 Leading order solution in the yielded regionTo the leading order (3.3) and (3.4) lead to0 =∂U0∂X+∂V0∂Y0 = −∂P0∂X+∂τXY,0∂Y0 = −∂P0∂YWe consider only Y > 0, as the flow is symmetric about the centreline.Integrating the X-momentum equation and satisfying (3.6a) gives∂U0∂Y+B sgn(∂U0∂Y) =∂P0∂X(Y − h) (3.9)The leading order position of the yield surface is denoted by Y = yy. Weknow at the interface of plug region (Y = yy) shear stress equates the yieldstress:∂P0∂X(Yy − h) = B sgn(∂U0∂Y) = −B. (3.10)273.3. O(δ) solution in the yielded layerIntegrating one more gives the velocity:U0(Y ) =up − up,0 if |Y | ≤ yy,up[1−(|Y | − yy)2(1− yy)2]− up,0 if |Y | > yy,(3.11)where the leading order plug velocity is given by:up =B2(yy − h(X))(1− yy)2. (3.12)In order to find yy, we impose the flow rate constraint across the channel,i.e. ∫ 1−1U0 dY = 2− 2up,0.Note that the leading order velocity in the droplet is given by continuity atthe interface as: U0 = up − up,0. Integrating across the channel gives us:y3y − yy(3 +6B) + 2 + 6hB= 0. (3.13)This cubic equation is similar to (2.16) and gives yy as a function of h(X)and B. It follows that up is also a function of both B and x. The dependenceof up on X via h(X) contradicts the idea that the plug region is rigid, i.e. thisis the so-called lubrication paradox of [42] discussed earlier in § 1.1.2. Theinconsistency is resolved at the next order in δ.3.3 O(δ) solution in the yielded layerWe focus here on the yielded layer of encapsulating fluid, avoiding the trickyissue of the first order momentum balance within the unyielded plug. Theidea is to correct the leading order solution, to account for the varying plugspeed predicted. In this, we shall assume that yy(X) is an O(δ) approxima-tion to the position of the true yield surface at Y = yT . The X-momentumbalance at 1st order in the yielded layer is:∂P1∂X=∂2U1∂Y 2(3.14)which requires two boundary conditions. One of these is no slip conditionat the wall (U1(Y = 1) = 0). To obtain the other condition notice that aty = yT we haveU(X, yT ) = up,0 = U0(X, yT ) + δU1(X, yT ) +O(δ2) (3.15)283.3. O(δ) solution in the yielded layerprovided that the spacing between droplets is long enough and the plugremains unyielded (recall that we have subtracted the pure fluid plug velocityup,0 from the X-component of velocity in translating to the moving frame).(3.15) can be expanded around yT and written asup,0 = U0(X, yy)+δ∂U0∂Y(X, yy)+δU1(X, yy)+δ2∂U1∂Y(X, yy)+O(δ2) (3.16)We can argue that the second term in the right hand side of the equation isat most of O(δ2), because:δ∂U0∂Y(X, yy) = δ∂U0∂Y(X, yT ) + δ2∂2U0∂Y 2(X, yT ) = δ2∂2U0∂Y 2(X, yT )Note that we have imposed the condition at yy(X), which is known, ratherthan at yT (X). However, assuming that yy(X) − yT (X) = O(δ), the dis-crepancy due to imposing at yy(X) is only at second order. Finally, thecondition on U1 is simply:U1(X, yy) =1δ(up,0 − up(X)) ≡ η. (3.17)This same value is adopted for all Y < yy in order to correct the plugvelocity. We will show in §A.3 that η = O(1).We integrate (3.14) impose the two boundary conditions and then imposethe flow rate constraint (∫ 1−1 U1 dY = 0) to find:∂P1∂X= −6ηyy + 1(yy − 1)3(3.18)andU1(X,Y ) = −η(Y − 1)(3Y + 3Y yy − 1− 4y2y − yy)(yy − 1)3if Y > yy (3.19)The last step is to find yT . We find yT by imposing zero shear rate at Y = yTand expanding about Y = yy with respect to δ:∂U0∂Y(X, yT )+δ∂U1∂Y(X, yT ) = 0 =⇒ yT (X) = yy(X)−δ∂U1∂Y (X, yy(X))∂2U0∂Y 2 (X, yy(X)).(3.20)An example of the perturbation solution is plotted in Figure 3.1a, show-ing U0(X,Y ) and yY (X) for an elliptical droplet with H = 0.2 and δ = 0.1293.4. Size of the yield surface perturbation from yy,0 and effects of inertiaFigure 3.1: Contours of the X-component of velocity obtained by the per-turbation solution; H = 0.2, δ = 0.1 and B = 1: a) leading order U0(X,Y );b) corrected velocity U0(X,Y ) + δU1(X,Y ). Uncorrected yield surface po-sition yy(X) in dashed line and corrected yield surface yT (X) in solid line.Pressure gradient for both solutions is also plotted.at B = 1. Although H and δ are relatively large here (i.e. for what mightbe considered asymptotic), the idea is simply to amplify visually the fea-tures of the perturbation solution. Figure 3.1b shows the corrected solutionU0(X,Y ) + δU1(X,Y ) and yT (X). The corresponding variations in pressuregradients are also shown on the top of each graph.3.4 Size of the yield surface perturbation fromyy,0 and effects of inertiaHaving now derived an expression for yT (X) we wish to examine the sizeof the yield surface perturbation from the Poiseuille flow yield surface, yy,0.As previously discussed in §2.5, it is this quantity that dictates the size ofthe inertial terms that we have neglected so far.We note that our perturbation solution has been derived under the as-sumption that h(X) ∼ O(δ) and we now define h¯(X) = h(X)/δ. Thecalculation proceeds in 2 steps. First we examine the departure of yy fromyy,0. We write yy = yy,0 + δξ1 and show that ξ1 is order 1. Substituting303.4. Size of the yield surface perturbation from yy,0 and effects of inertiaboth yy and h¯ into (3.13) and expanding in powers of δ:O(1) : y3y,0 − yy,0(3 + 6/B) + 2 = 0, (3.21a)O(δ) : 3y2y,0ξ1 − ξ1(3 + 6/B) + 6h¯B= 0, =⇒ ξ1 =h¯1 + 0.5B(1− y2y,0)(3.21b)The first expression is simply the Buckingham equation (2.16), as ex-pected. The second expression defines ξ1. Now, since up is defined byyy and h, straightforward expansion of (3.12) in powers of δ verifies thatη = O(1); see Appendix A.3.Secondly, combining(3.20) and (3.21b)yT (X)− yy,0 = δ(ξ1(X)−∂U1∂Y (X, yy)∂2U0∂Y 2 (X, yy)). (3.22)and show that (algebraic details can be found in Appendix A.3):yT − yy,0 = δh¯(B(yy,0 + 1)(y3y,0 − yy,0(3 +6B ) + 2)yy,0(yy,0 − 1)(By2y,0 −B − 2))+O(δ2). (3.23)We see that the leading order term on the right-hand side vanishes dueto (3.21a), leading to: yT − yy,0 ∼ O(δ2). A similar exercise shows thatU0 + δU1 = Up(y) + O(δ2), i.e. the velocity field is given to O(δ2) by thePoiseuille solution of the droplet free flow (algebraic details can be found inAppendix A.3).To summarise, insertion of an O(δ) droplet within the plug region causesonly anO(δ2) perturbation in both the yield surface position and the velocityfield. However, the stress perturbation is of O(δ), as is the pressure gradientperturbation. Additionally, the stress field within the unyielded plug regionremains indeterminate. To explore this surprising feature, Figure 3.2 plotsyT (X) and yy(X) for a range of different droplets at B = 10. As we seethe corrected yield surface position is nearly constant and coincides withyy,0 (equal to the end values of yy(X)). Considering now the question ofinertial effects, we see that for an O(δ) encapsulated droplet as considered,the uniform Poiseuille flow gives an O(δ2) approximation to the velocity fieldin the region of the droplet (as well as away from the droplet). The yieldsurface perturbation is also of O(δ2). This suggest that the appropriatevariation in the streamlines away from the uniform Poiseuille solution in thepresence of an encapsulated droplet is in fact O(δ2) (instead of the O(δ)used to derive and construct our perturbation approximation), and thatconsequently inertial effects may be legitimately neglected for δ2Re 1.313.4. Size of the yield surface perturbation from yy,0 and effects of inertia−0.05 −0.025 0 0.025 0.0500.−0.1 −0.05 0 0.05−0.2 −0.1 0 0.1−0.05 −0.025 0 0.025 0.0500.−0.1 −0.05 0 0.05−0.2 −0.1 0 0.1−0.05 −0.025 0 0.025 0.0500.−0.1 −0.05 0 0.05−0.2 −0.1 0 0.1 b b(a) (b) (c)(d) (e) (f)Figure 3.2: Position of the yield surfaces yT (X) (solid line) and yy(X)(dashed line) for a range of different (elliptical) droplets at B = 10. Toprows L = 0.5, bottom row L = 0.75 and columns left to right H = 0.05, 0.1and 0.2. Note that yy,0 is equal to the values of yy(X) at two ends. Due tothe adopted scaling, all elliptic droplets are mapped to a circle with radiusof H.32Chapter 4Numerical SimulationIn this chapter we describe the detail of numerical analysis that we have con-ducted to solve the encapsulation problem without having any geometricalconstrain on the drop’s shape.4.1 Assumptions and simplifications.Motivated by the asymptotic results, we proceed under the assumption thatan encapsulated droplet held within the plug will not significantly perturbthe yield surface, and consequently that the inertial terms may be smalleven for significant Re. This we verify a posteriori from our numericalresults. We also consider that the droplet and encapsulating fluid have thesame density, i.e. χ = 0.As we consider a Stokes flow problem, we impose symmetry conditionsat x = 0, x = −l/2 and y = 0, and solve only for: (x, y) ∈ [−l/2, 0]× [0, 1];(see Figure 4.1). In the encapsulating fluid the Stokes equations are0 =∂u∂x+∂v∂y, (4.1a)0 = −∂p∂x+∂τxx∂x+∂τxy∂y, (4.1b)0 = −∂p∂y+∂τxy∂x+∂τyy∂y. (4.1c)The constitutive laws are (2.10) & (2.12). The droplet domain is y ≤ h(x)and within the droplet the fluid is assumed to move uniformly at the speedof the unyielded plug. Thus, we solve only in the encapsulating fluid domain;334.1. Assumptions and simplifications.Figure 4.1: Schematic description of the problemsee Figure 4.1. The boundary conditions are:− p(0, y) + τxx(0, y) = 0, v(0, y) = 0, y ∈ [H, 1] (4.2a)u(x, 1) = 0, v(x, 1) = 0, x ∈ [−l/2, 0] (4.2b)− p(−l/2, y) + τxx(−l/2, y) = Gl/2, v(−l/2, y) = 0, y ∈ [0, 1] (4.2c)τxy(x, 0) = 0, v(x, 0) = 0, x ∈ [−l/2,−L] (4.2d)τnt(x, h(x)) = 0, − p(x, h(x)) + τnn(x, h(x)) = 0, x ∈ [−L, 0] (4.2e)Note in (4.2b) that for simplicity we compute in a fixed frame of referenceinstead of the moving frame, as translating by a constant up,0 can be im-posed afterwards with no loss of generality. By imposing (4.2c) instead ofthe uniform Poiseuille flow we allow for the fact that the encapsulating fluidbetween droplets may not be fully developed. Indeed, below we study therequired droplet spacing parameter l, needed for the droplets not to interactwith each other. The parameter Gl/2 in (4.2c) is essentially the pressuredrop along the computational domain. Each computation involves an iter-ation on G to ensure that the mean velocity is unity, i.e. we satisfy∫ 10u(x, y) dy = 1, (4.3)which has a unique solution due to the monotonicity of the flow rate withrespect to the pressure drop.344.2. Rheolef4.2 RheolefAs mentioned earlier, regularisation methods do not predict the position ofthe yield surfaces with certainty. Since we investigate whether or not thedroplet is fully encapsulated in unyielded fluid, it is essential to use aug-mented Lagrangian method to distinguish the true plug from regions withmerely small strain rate. This necessity eliminates the options of using acommercial CFD software. The augmented Lagrangian algorithm is imple-mented within the code Rheolef, a C++ Finite Element Method library byP. Saramito and colleagues in Grenoble.6 In order to illustrate how powerfulRheolef is, here we show the first example of its tutorial. This emphasisesthat working with Rheolef eliminates all the hassle of code compiling. Thefollowing code solves 4u = −1 in Ω with u = 0 on ∂Ω when Ω is a givencomputational domain.Program 4.1 Solving Poisson’s equation with Dirichlet boundary conditionin Rheolefint main (int argc, char** argv) {environment rheolef (argc, argv);geo omega (argv[1]);space Xh (omega, argv[2]);Xh.block ("boundary");trial u (Xh); test v (Xh);form a = integrate (dot(grad(u),grad(v)));eld lh = integrate (v);eld uh (Xh);uh ["boundary"] = 0;solver sa (a.uu());uh.u = sa.solve (lh.u());dout << uh;}4.3 Mesh adaptationA particular feature of our implementation in Rheolef is intensive meshadaptivity, which focuses the mesh around the yield surfaces. It has been6for more details visit http://www-ljk.imag.fr/membres/Pierre.Saramito/rheolef/354.4. Convergenceshown by [67] that the rate of dissipation defined byϕ = ||γ˙||2 +B||γ˙|| (4.4)exhibits a sharp jump across the yield surface. This behaviour can be utilizedin order to find the approximate location of plug interface and produce afiner mesh nearby. Details of implementation can be found in [67, 64, 61, 65].Starting with an initial mesh, with a typical mesh size of say d0 ≈ 0.02, wecompute the velocity and pressure using AL method. Then we calculate ϕand regenerate the mesh such that the gradient of ϕ is equally distributedin all cells. This process is repeated typically four times until we get adesirable mesh configuration. Figure 4.2 shows cycles of the adaptation forone specific case. After the fourth cycle the typical mesh sizes are: in theplug: 0.02; near the yield surface: 0.001× 0.004; near the wall: ConvergenceStarting with a smaller d0 results in a progressively smaller mesh at eachstep of the adaptivity and hence a better converged solution. The effects ofadaptivity on convergence are also significant for the first few cycles. Thisis illustrated in Figure 4.3 where we have compared the numerical resultsof the code for plane Poiseuille flow with the exact solution from §2.4. Thisshows that ||eu|| and ||ep||, defined below, decay rapidly over the first fewcycles of adaptation before leveling out, i.e. reducing the error after 4 cyclesrequires a smaller initial mesh d0.||eu|| =1Ncells√∑Ncells|ucomp − uexact|2 (4.5a)||ep|| =1Ncells√∑Ncells(pcomp − pexact)2 (4.5b)We are also able to compare our computed results against the asymptoticsolutions over a range of small δ, e.g. fixed small H and varying L. Thistype of comparison verifies that the errors of leading order and correctedfirst order solutions are O(δ) and O(δ2), respectively. The range of δ testedis however limited since for very small H the droplet becomes a cut inthe encapsulating fluid, increasingly singular at the ends. On the otherhand for larger L the encapsulating plug region breaks. As an aside, we364.4. Convergence−1 −0.75 −0.5 −0.25 bc d−0.35 −0.25 −−0.27 −0.260.710.720.73Figure 4.2: Mesh adaptation cycles for the case of B = 20, H = 0.1 andL = 0.5, where we start with an initial mesh scale, d0 ≈ 0.02: a) cycle 1,3165 elements; b) cycle 2, 12689 elements; c) cycle 3, 21425 elements; d)cycle 4, 37290 elements. After the 4th cycle the typical mesh sizes are: inthe plug: 0.02; near the yield surface: 0.001× 0.004; near the wall: 0.002.374.4. Convergence0 0.5 1 1.5 2x 10510−910−810−710−610−510−410−3Number of Cells||e u||  cycle 1Adaptation cycle 0(a)cycle 3 cycle 4cycle 20 0.5 1 1.5 2x 10510−510−410−310−210−1100Number of Cells||e p||  (b)cycle 4cycle 3cycle 0cycle 1cycle 2Figure 4.3: Code Validation: Error of numerical results for channel Poiseuilleflow compared to the exact solution; © : B = 10, ♦ : B = 20 and 4 : B =50. a) velocity and b) pressure384.5. Effects of droplet spacingFigure 4.4: Example of the effects of droplet spacing for B = 20, H = 0.5and L = 0.6: a) l = 1.625; b) l = 1.8; c) l = 2.5. The colour map shows therange of speeds in the flow and grey areas show unyielded fluid.mention that the plane Poiseuille flow solution also gives an O(δ2) errorwhen compared against the computed results and interestingly a smallererror than the corrected perturbation method.4.5 Effects of droplet spacingAs we have already seen, the length of droplets has a significant effect onbreaking of the encapsulating plug. With a restriction in L, if one is in-terested in the throughput of encapsulated droplets it is necessary to un-derstand the effect of droplet spacing on the breaking process. Figure 4.4shows an example sequence of computations for which the separation dis-tance l− 2L drops below a critical value and the plug regions start to yield.In this case, for l = 1.625 (Figure 4.4a) the plug has broken whereas forlarger l an intact region exists. If the plug yields around the droplet thecomputed solution remains valid as solution to the Stokes flow, but in termsof an evolution problem the interface is free to deform and we would expectthis deformation to induce a velocity field within the droplet.Apart from the issue of plug breaking, a separate consideration is whetherl is long enough for the flow to become fully developed between the droplets.394.5. Effects of droplet spacing1 2 3 4lmax(2H,2L)1 2 3 41015202530lmax(2H,2L)∆p ∆x  (a) (b)Figure 4.5: Variation in the pressure drop upstream of the droplet plottedversus a scaled spacing between droplets. a) B = 10, O: H = 0.45, L = 0.6;©: H = 0.55, L = 0.15; : H = 0.05, L = 0.75; ?: H = 0.25, L = 0.25and ♦: H = 0.375, L = 0.5 b)B = 20, O : H = 0.5, L = 0.6; ©: H = 0.65,L = 0.2; : H = 0.1, L = 0.9; ?: H = 0.25, L = 0.25 and ♦: H = 0.375,L = 0.5. The broken lines show the pressure drop for Poiseuille flow.The method we have chosen to study this compares the mean pressure gra-dient behind the drop with the analytical solution for plane Poiseuille flow:∆p =pl − pLl − L,where pl and pL are the pressures at x = −l/2 and x = −L, respectively;see Figure 4.5. For a range of droplet shapes we observe that at sufficientlylarge l the encapsulation of the droplet is intact and the ∆p approachesthe analytic value for a yield stress fluid. In fact ∆p is always very slightlybelow the analytical value, due to a reduction very close to x = −L, wherepresumably the fluid senses the presence of the droplet. As a rule of thumb,if l is at least 2.5 times max{H,L} we find that the pressure drop valuesremains near constant.40Chapter 5Results and Discussion5.1 Encapsulation and failureWe have computed approximately 400 encapsulation flows in channel ge-ometry, as described. These cover the approximate ranges: 5 < B < 200,0.025 < H < 0.9 and 0.025 < L < 1.25. Figure 5.1 illustrates some ofthe general trends observed. For each row of computations the Binghamnumber is held constant. The left-hand column has L kept constant and Hvaried; the right column has variable L as H stays fixed. The colourmapplots the speed of the fluid, with the grey area denoting the unyielded plug.We observe that there essentially are two different regimes of plug breaking:Slender droplets : If the droplet aspect ratio δ = H/L is small (in otherwords, the length of drop is much larger than its width) then theplug appears to yield initially close to the two tips of the droplet (leftcolumn of Figure 5.1).Fat droplets : If the width of the droplet is similar or larger than thelength, the plug appears to yield closer to x = 0, in the widest part ofthe droplet (right column of Figure 5.1).Essentially this suggests that increasing either L or H will eventually yieldthe plug around the droplet, as might be expected.The other remarkable observation is that the yield surface position isapproximately constant for the simulations shown in Figure 5.1, up to anddirectly after the plug breaks (although of course with breaks in it afteryielding). This implies that the stress distribution to the outside of theplug is not greatly affected by plug breaking. Thus, the deduction fromthe asymptotic analysis of the past section: of an O(δ2) perturbation of theyield surface due to the droplet, appears to be more widely valid even forlarger aspect ratios δ = H/L.To examine the breaking mechanism, we have plotted the various stressesin Figure 5.2 for one sequence of slender drops (B = 20, H = 0.2 and in-creasing L). The first noticeable feature is that as L increases through the415.1. Encapsulation and failure(a)Figure 5.1: Speed distribution (u) as the droplet gets larger. Rows representconstant Bingham number; H = 0.2 over left column (slender drop) andL = 0.5 over right column (fat drop). a) B = 5 and L = 0.55, 0.6, 0.65,0.675 and 0.7; b)B = 5 and H = 0.3, 0.34, 0.375, 0.4 and 0.425; c)B = 20and L = 0.8, 0.85, 0.9, 0.925 and 1;d) B = 20 and H = 0.5, 0.55, 0.575, 0.61and 0.625; e)B = 50 and L = 0.95, 1, 1.05, 1.15, 1.2 f)B = 50 and H = 0.6,0.65, 0.69, 0.725 and 0.75 425.2. Maximum size of encapsulated dropletsbreaking transition the distribution of stresses remains qualitatively similar.Thus, the breaking itself does not affect the overall pattern of stresses whichare dominated primarily by changes in droplet shape. Upstream and down-stream of the droplet the stresses develop relatively quickly (within a chan-nel width) and adopt distributions close to fully developed. The verticallyoriented yielded regions (cracks?) coincide with a region within which theshear stress is focused: τxy becomes large and negative near the tips of thedroplet. This is combined with significant extensional stresses τxx = −τyy.The basic mechanism appears to be that outlined in §2.6. In the absenceof extensional stresses, the pressure must transit from its shear-layer valuesto zero at the droplet surface. This driving force is evidently larger at thetips of the droplet. One way of reducing the burden on the change in pis by taking up a part of the pressure within the extensional stresses. Forexample, on the interface where it is nearly flat the normal stress conditionis p = τyy and we can observe that |τyy| becomes significant, changing signtowards the yield surface. In the centre of the droplets it appears that thedistribution of τxy is shifted upwards by h(x), but near to the ends of thedroplet the stresses adapt to the far-field in a relatively narrow region.The stress distributions in the case of fat droplets are illustrated in Fig-ure 5.3. for B = 20, L = 1.4 and increasing H. Although the far-fieldsare similar to the slender droplet, the mechanism of breaking appears to bedifferent. Again we have a transition as the geometry is increased. How-ever, this time as H increases, note that from the momentum balance (andobservation that the yield surface position is unchanged), we expect∂τxy∂y∼ O(Byy −H),in the centre of the droplet (within a layer that is getting progressivelyshort). If we suppose that the pressure gradient is partly constrained by thefar-field plane Poiseuille flow pressure gradient, then we might expect thatthe y-derivative of τxy is balance by x-derivatives in τxx which we observe inFigure 5.3. The yielding appears to occur at the ends of the central shearlayer.5.2 Maximum size of encapsulated dropletsBy successively iterating with respect to H and L we are able to approxi-mately determine the limiting size of droplet that can remain encapsulatedwithin the plug region. The procedure that we took here was to simulate435.2. Maximum size of encapsulated dropletsFigure 5.2: Stress distribution as the droplet gets larger. B = 20, H = 0.2and; L = 0.8, 0.85, 0.9, 0.925 and 1. a) speed; b)pressure; c)τxy; d)τxx;e)σxx = −p+ τxx and f)σyy = −p+ τyy445.2. Maximum size of encapsulated droplets(b)(f)Figure 5.3: Stress distribution as the droplet gets larger. B = 20, L = 0.7and H = .5, 0.55, 0.575, 0.61 and 0.625. a) speed; b)pressure; c)τxy; d)τxx;e)σxx and f)σyy455.2. Maximum size of encapsulated droplets0 0.2 0.4 0.6 0.8 1 1.2 1.400.  Yielded PlugUnyielded PlugFigure 5.4: Maximum size of encapsulated droplets, computed for five dif-ferent B. ©: B = 5; : B = 10; 4: B = 20, ♦: B = 50 and ?: B = 200flow around a drop and find the plug region. If the plug region is intact, weenlarge the drop until the plug starts to yield. The accuracy of incrementsis 0.01 in both directions, H and L. Although this process might seemstraightforward, some issues have to be dealt with. To explain these issues,note that to determine the surface of plug regions, we could take two similarapproaches: calculate either strain rate (γ˙) or stress (τ). In either case acomparison with a constant value (0 or B, respectively) is necessary here toidentify the surface of plug regions. Generally comparing equality is not anice operation in computational analysis, due to small numerical errors. Forexample, assume we have a computational geometry with typical mesh sizeof d. If the round-off error is 1, then the order of error we might have forthe strain rate is 1d . Thus, we need to allow for such errors in finding theactual surface of plug region and our comparison algorithm is as follows:{|γ˙| < 2 or |τ | < B(1 + 2) plug regionotherwise yielded regionwhere 2 =1d .Figure 5.4 shows the critical values in the (L,H)-plane for five different465.3. Encapsulation with fluids of different densitiesvalues of Bingham number, B. For a given Bingham number, any drop lyingbelow the associated curve does not break the plug region. We see that forany Bingham number (no matter how large) we have a maximum value forthe length and height of the droplet. The maximum height is mostly limitedby the plug width. Although there is no similar simple physical limit for themaximum length, we still see that the growth of stresses for long dropletsyields the plug. In fact the shape of the curves at small H indicates thatfor slender droplets the encapsulation becomes very sensitive to changes inthe length, i.e. small increases in the length can break the plug (see alsoFigure 5.2a, c & e).A last observation here is that the increase in size of droplet with Bis much less than linear. Therefore, just increasing the yield stress of theencapsulating fluid does not allow correspondingly large droplets. The un-derlying reasons why seem to be captured in the simple analysis of §2.6.Long droplets break near the ends. Breaking arises from the need to bridgethe discrepancy between non-zero (frictional) pressures in the yielded fluidlayers and zero normal stress on the droplet surface. This discrepancy in-creases with L but also with B, i.e. the frictional pressure scales like xB/yy,0.5.3 Encapsulation with fluids of differentdensitiesWe now explore the effect of allowing a density difference between the fluids.Assuming the encapsulated fluid remains in steady rigid motion, the pressurefield with in the droplet is: p = χx+ constant, where we recall that χ =(ρˆd− ρˆ)gˆRˆ2/(µˆUˆ0) represents the competition between buoyant and viscousstresses. Considering the preceding discussion and the simple analysis in§2.6, an intuitive notion is that by selecting χ = −B/yy,0 the pressuregradient in the droplet would match that in the yielded fluid layers. Thiswould eliminate the need for normal stress gradients across the plug, andhence presumably reduce τ .We test the above notion computationally. The code described earlieris adapted by discretizing both the droplet and carrier fluid domains. Theadditional body force term χ can be added to the droplet x-momentumequation, (or instead −χ added to the carrier fluid equations). Differentfluid properties are assigned in each domain (e.g. B = 0 is inside the droplet)and the Stokes flow solution is computed on both domains simultaneouslyusing the augmented Lagrangian algorithm described earlier.Firstly, we address the question of whether the density difference might475.3. Encapsulation with fluids of different densitiessignificantly affect breaking of the plug region. Figure 5.5a shows the max-imum height H before the plug yields, for two droplets with different Lat B = 20. We observe that the shorter droplet (L = 0.4) is relativelyinsensitive to χ, whereas the longer droplet (L = 1.75) shows remarkablesensitivity. Sensitivity to L per se is to be expected as the effect of χnaturally amplifies with x (or L), due to the channel orientation. We seethat if the droplet is heavier than the carrier fluid (χ > 0), encapsulationfails at a smaller sizes for the drop. On the other hand, for lighter droplet(χ < 0) buoyancy tends to stabilize the encapsulation, allowing larger dropsto be successfully encapsulated. For the case shown in Figure 5.5a, in whichB = 20, we see that maximum size of the drop happens at χ ≈ −27 to −28.This approximately coincides with the frictional pressure gradient of theplane Poiseuille flow, (for B = 20 we find ∂p∂x = −Byy,0≈ −27.8). Thus, thenotion that selecting χ ≈ −B/yy,0 may provide an optimal size of encapsu-lated droplet appears reasonable. Essentially buoyancy acts to balance thefrictional pressure losses for χ > −B/yy,0, but on decreasing χ still furtherstrong buoyancy forces act to exert stress on the plug.Note that although the normal stress discrepancy may be reduced forχ ≈ −B/yy,0, the tangential stresses in the plug are still perturbed by thepresence of the droplet. Thus, the droplet may still break the plug. Figure5.5b shows the critical values of H and L for B = 20, above which the plugyields for three different χ = −20 , 0 and 20. Again we see favourable densitydifference (lighter droplets) allows the encapsulation of longer drops, whereasunfavourable density difference (heavier droplets) weakens the encapsulationcapabilities.Figure 5.6 explores variations in the stress fields as the density of thedroplet is varied. The plug close to the central part of the droplet is largelyunaffected by χ, but on increasing |x| we observe generation of extensionalstresses towards the ends of the droplet within the plug (Figure 5.6d). In-creasing χ amplifies the size of extensional stresses and of the pressure inthese regions, (Figure 5.6b & d). The stresses combine so that σyy showslittle variation in y (Figure 5.6f), which implies that the pressure variationin y in negatively amplified in τxx. The net effect is that the front of thedroplet is stretched and the rear is compressed (in the x-direction); see Fig-ure 5.6e. These normal stress gradient results in significant shear stresses(Figure 5.6c). The underlying pattern is not changed by χ > 0, but simplyamplified.485.3. Encapsulation with fluids of different densities−80 −60 −40 −20 0 20 4000.χH  0 0.5 1 1.5 2 2.500.  Yielded PlugUnyielded PlugFigure 5.5: B = 20 (a) Variation of maximum height (maximum H whichdoes not break the plug) for two different length of drop : L = 0.4 and©: L = 1.75 with density difference (χ). (b) Maximum size of drop for forthree different χ. ©: χ = 20; : χ = 0 and 4: χ = −20.495.3. Encapsulation with fluids of different densitiesFigure 5.6: Stress distribution as the droplet becomes heavier: B = 20,L = 0.9 and H = .2. From top to bottom in each panel: χ = −10, −1, 0, 1and 10. a) speed; b) pressure; c) τxy; d) τxx; e) σxx and f) σyy.50Chapter 6Axisymmetric GeometryFor the purposes of practical application we now consider the analogousflow to that described in §2.1, but within a vertical pipe of radius Rˆ. Wedescribe the problem in a cylindrical coordinate system, assuming that theflow is axisymmetric with droplets traveling along the pipe axis. As in §2.1 amoving coordinate system is adopted. Lengths and velocities are scaled withRˆ and Wˆ0 respectively where Wˆ0 is the mean velocity along the pipe, i.e. theflow rate is piRˆ2Wˆ0. Denoting by Wˆp the fully developed plug velocity ofthe Hagen-Poiseuille flow (with no droplets), the variables are as follows:(r, z, θ) =(rˆRˆ,zˆ − WˆptˆRˆ, θ); (u,w) =(uˆrWˆ0,uˆz − WˆpWˆ0), (6.1)and t = Wˆptˆ/Rˆ. The deviatoric stresses are scaled with µˆWˆ0/Rˆ as is themodified pressure (on subtracting the static pressure of the carrier fluid).pˆ = pˆ0 + ρˆgˆzˆ +µˆWˆ0Rˆp (6.2)Exploiting axisymmetry of the geometry, the dimensionless form of con-tinuity and momentum equations governing carrier fluid are1r∂∂r(ru) +∂w∂z= 0 (6.3a)Re(u∂w∂r+ w∂w∂z)= −∂p∂z+1r∂∂r(rτrz) +∂τzz∂z(6.3b)Re(u∂u∂r+ w∂u∂z)= −∂p∂r+1r∂∂r(rτrr) +∂τrz∂z−τθθr(6.3c)The constitutive relations are similar to (2.10), but now γ˙ =√γ˙2rr + γ˙2zz + γ˙2θθ + 2γ˙2rzis the rate of strain, and τ =√τ2rr + τ2zz + τ2θθ + 2τ2rz is the deviatoric stress.The components are given by:γ˙rr = 2∂u∂r; γ˙θθ = 2wr; γ˙zz = 2∂w∂z; γ˙rz = γ˙zr =∂w∂r+∂u∂z. (6.4)516.1. Poiseuille flow solutionReynolds number and Bingham number are defined as Re = ρˆWˆ0Rˆ/µˆ andB = τˆY Rˆ/µˆWˆ0.As before, we consider a single droplet from our periodic train of droplets.The scaled droplet length is L and the maximal radius is H. The dropletsurface is steady within the moving frame of reference and is denoted r =h(z). Within the droplet we must satisfy:1r∂∂r(ru) +∂w∂z= 0 (6.5a)φRe(u∂w∂r+ w∂w∂z)= −∂p∂z+m1r∂∂r(r∂w∂r) +m∂2w∂z2+ χ (6.5b)φRe(u∂u∂r+ w∂u∂z)= −∂p∂r+m1r∂∂r(r∂u∂r) +m∂2u∂z2−mur2(6.5c)where χ = (ρˆd−ρˆ)gˆRˆ2Wˆ0µˆ. Similar boundary and interfacial conditions are satis-fied as for the plane channel problem earlier. As before, it is assumed thatthe droplet translates uniformly within the moving plug.6.1 Poiseuille flow solutionTaking a similar approach to that of §A.1, we can readily find that in a fixedframe of reference the Hagen-Poiseuille solution in pipe geometry is resolvedby:Wp(r) =wp,0 if r ≤ ry,0,wp,0[1−(r−ry,01−ry,0))2]if r ≥ ry,0,(6.6)where wp,0 = B2ry,0 (1−ry,0)2 is the dimensionless plug velocity. On imposingthe flow rate constraint we find that the yield surface position satisfies thefollowing Buckingham equation:r4y,0 − 4ry,0(1 +3B) + 3 = 0, (6.7)which has a unique zero in the interval (0, 1).6.2 Slender dropFor practical reasons it is advantageous to know if encapsulated dropletswithin a pipe flow also result in vanishing small perturbations of the flow526.2. Slender dropstreamlines. Thus, similar to §3.1 we re-scale the problem assuming a smallslender droplet:r = R, θ = Θ, zδ = Z, u = δU, w = W, pδ = P. (6.8)Upon rescaling and neglecting the inertial terms, the governing equationsare:1R∂∂R(RU) +∂W∂Z= 0 (6.9a)0 = −∂P∂Z+1R∂∂R(RτRZ) + δ2∂τZZ∂Z(6.9b)0 = −∂P∂R+ δ21R∂∂R(RτRR) + δ2∂τRZ∂Z− δ2τΘΘR(6.9c)where τrz = τRZ , τrr = δτRR, τθθ = δτΘΘ and τzz = δτZZ . The constitutiverelations are similar to (2.10), except that now:τ =√τ2RZ + δ2τ2RR + δ2τ2ΘΘ + δ2τ2ZZ (6.10a)γ˙ =√γ˙2RZ + δ2γ˙2RR + δ2γ˙2ΘΘ + δ2γ˙2ZZ (6.10b)andγ˙RZ =∂W∂R+ δ2∂U∂Z(6.11a)γ˙RR = 2∂U∂R(6.11b)γ˙ΘΘ = 2WR(6.11c)γ˙ZZ = 2∂W∂Z(6.11d)We adopt a regular perturbation expansion in δ: (P,W,U) = (P0,W0, U0)+δ(P1,W1, U1) + · · · , substitute into the equations and derive the followingleading order problem:1R∂∂R(RU) +∂W∂Z= 0 (6.12a)0 = −∂P∂Z+1R∂∂R(RτRZ) (6.12b)0 = −∂P∂R(6.12c)536.2. Slender dropIntegrating the axial momentum equation from R = h yieldsτRZ =∂W0∂R+Bsgn(∂W0∂R)=∂W0∂R−B =12∂P∂Z.(R−h2R)(6.13)and defining the yield surface R = Ry as where τRZ = −B we find:∂W0∂R=12∂P∂Z(R−Ry −h2R+h2Ry):∂P∂Z= −2BRyR2y − h2(6.14)We can reconstruct the velocity in the fixed frame, using the no-slip condi-tion at R = 1:W0(r)+wp,0 ={wp if 0 ≤ R ≤ Rywp[1−(0.5(R−Ry)2+h2(R−Ry)/Ry+h2 ln(Ry/R)0.5(1−Ry)2+h2(1−Ry)/Ry+h2 lnRy)]if Ry ≤ R ≤ 1(6.15)where the plug velocity is given bywp = BRyR2y − h2((1−Ry)22+h2Ry(1−Ry) + h2 lnRy)(6.16)Finally we use the unit flow rate condition to find the pressure gradient(equivalently Ry). This leads to the following Buckingham equation:R4y − 4(1 + 3/B)Ry + 3 +2h2Ry[(1−Ry)2(2 +Ry) +6B]= 0, (6.17)which may be solved numerically for Ry.Comparing the solution above with that in §6.1, we note that the lead-ing order plug velocity, velocity profile and the Buckingham equation areall O(h2) perturbations from the unperturbed Poiseuille flow solution, c.f.(6.15) , (6.17), (6.6) & (6.7).The implications of this are two-fold. Firstly, in terms of the perturba-tion procedure there is no need to seek a correction to the solution at firstorder in δ. In practical terms, this means that for any small slender dropof aspect ratio δ the streamlines in the yielded layer will be perturbed byO(δ2), suggesting that inertial effects on the flow are of O(δ2Re) as for thechannel flow. Secondly, this suggests that a perturbation approximation ofthe above form would only require a first order correction when h ∼ δ1/2.As δ = H/L is the droplet aspect ratio, this suggests L ∼ δ−1/2. This hintsthat longer and wider droplets may be accommodated in the pipe geometry.546.3. Examples of encapsulation and failure6.3 Examples of encapsulation and failureWe now present computed solutions of droplets as either L orH are increaseduntil breaking point. Figure 6.1 shows yielding of the “slender” droplet andFigure 6.2 shows the “fat” droplet. Superficially, the yielding mechanismare similar to the channel flow in where the plugs break (ends or near thecentre). However, as suggested by the previous section, the plugs generallybreak at larger values of L.The main reason is in the normal stresses. The pressure distributionsfound are qualitatively similar between pipe and channel geometries. At theyield surface in the channel −σyy is continuous and approximately equal tothe pressure. The normal stress must then vanish at the droplet surfacewhich transfers some of the pressure into τyy. Since τxx + τyy = 0 we ex-perience extensional stresses along the plug, i.e. τxx 6= 0, and particularlynear the ends of the droplet. These normal stress gradients promote shearstresses τxy, but additional shear stresses result directly in the widest partof the droplet where the stress gradients must increase.The pipe differs firstly in that although normal stresses are generated atthe yield surface, τrr can now be balanced by the hoop stress τθθ as well asby τzz. Indeed, comparing Figs. 6.1e & f, we see that τrr and τθθ largelycancel out over much of the plug, leaving only small bands where τzz issignificant (Figure 6.1d). As a result, the shear stresses in much of the plugappear drastically reduced; see Figure 6.1c, occurring mostly close to z = 0.6.4 Maximum size of encapsulated dropletsFollowing approximately 300 computations we have been able to quantify themaximal size of encapsulated droplets at four different Bingham numbers;see Figure 6.3. In order to have a better comparison between channel andpipe geometries, we include the channel data for B = 10. We see that inthe range of slender droplets encapsulation is much easier in pipe geometrythan in the channel although for short droplets the channel allows larger H.The computed results for the slender droplets are suggestive of L → ∞ asH → 0. Although consistent with the notion that H ∼ δ1/2 and L ∼ δ−1/2(see §6.2), the slender geometries are not ideal for computations.556.4. Maximum size of encapsulated dropletsFigure 6.1: Stress distribution in the axisymmetric geometry as the dropletgets longer: B = 20, H = 0.4 and L = 1.2, 1.25, 1.3, 1.34 and 1.45. a)speed; b) pressure; c) τrz; d) τzz; e) τrr and f) τθθ.566.4. Maximum size of encapsulated droplets(c)Figure 6.2: Stress distribution in the axisymmetric geometry as the dropletheight is increased: B = 20, L = 0.5 and ;H = 0.4, 0.45, 0.5, 0.525 and 0.55.a) speed; b) pressure; c) τrz; d) τzz; e) τrr and f) τθθ.576.4. Maximum size of encapsulated droplets0 0.5 1 1.5 2 2.5 3 3.500.  UnyieldedPlugYielded PlugFigure 6.3: Maximum size of encapsulated droplets for four different Bing-ham numbers. ©: B = 5; : B = 10; 4: B = 20, ♦: B = 50 Includedfor comparison is a single curve showing the maximum size of droplet in thechannel geometry for B = 10 (broken line).58Chapter 7ConclusionsIn visco-plastically lubricated flows, multi-layer flows are kept stable throughthe yield stress of the lubricating fluid, which prevents the interface fromdeforming; see e.g. [21, 51, 38, 36, 37]. In [33] this concept has been extendedto multi-layer flows in which the interface is shaped, e.g. wavy, but againthe yield stress prevents interfacial deformation; see [33].In this thesis we have extended the concept in a different direction.Steady Poiseuille flows of yield stress fluids along uniform ducts are charac-terized by unyielded central plug regions, traveling at constant speed alongthe duct. Here we have investigated whether we may encapsulate dropletsof a second fluid within the plug region of the yield stress fluid, with againthe shape held constant by the yield stress of the fluid. The key advan-tage of the proposed methodology is that the yield stress holds the interfacerigid, instead of relying on capillary forces as is common in many dropletencapsulation studies, allowing Macro-size drops to be encapsulated.7.1 Summary of the resultsIn general, we have shown that the above methodology is feasible; i.e.Macro-size drops can be encapsulated inside the plug region of a poiseuilleflow of a Bingham fluid in both Cartesian (channel) and cylindrical (pipe)coordinates. In particular,• Introducing a droplet inside the plug region, travelling with the steadyspeed of plug does not deform the interface but does induce additionalstresses within the carrier fluid. It is the latter that may lead tobreaking of the plug.• In both geometries (channel and pipe) the yield surface of the plugregion containing the droplet is not significantly perturbed, apparentlyremaining near-planar as the droplet size is increased up to the pointat which the plug breaks. We observed this in both analytical andcomputational results.597.1. Summary of the resultsFigure 7.1: Encapsulation of drops with exotic shapes.• From the asymptotic solution it appears that the deformation of theplug surface is only of second order (O(δ2)) for slender droplets, whereδ = H/L 1 represents the droplet aspect ratio.• We can expect that inertial effects in the yielded part of the flow willscale like εRe, where ε is the aspect ratio of the streamlines, (i.e. ε = δ2for slender droplets). Thus, we expect that the solutions and resultswe have computed here will retain their validity as approximations tothe solution for significant Re, even though they are computed underStokes flow assumptions.• We know inertial terms are responsible for flow instability. There-fore, since our results suggest that the velocity perturbation due tothe droplet is small, there is the possibility that these flows may beobservable up to high Re.• As for failure of encapsulation, two main modes are uncovered: dropletsthat become increasingly slender tend to break the plug at the endsof the droplet, whereas those that become increasingly wide tend tobreak near the centre of the droplet. Loosely speaking, the end mode ofbreakage appears to be driven by normal stress development whereasthe central mode is driven by tangential stress gradients.607.1. Summary of the results• For elliptic (ellipsoidal) droplets we have computed the critical H &L at which breaking occurs for both channel and pipe geometries.Increasing the yield stress (via B) does increase the limiting dropletsize, but not as much as might be expected. This is because the yieldstress also increases the size of stresses found in the flow, i.e. consideredat constant flow rate.• The limiting values of H & L are different for channel and pipe, asmight be expected. At the same Bingham number for short drops thechannel may allow larger H for the same L. Here the limiting factoris anyway the yield surface position, which is similar as a function ofB. The limitation in terms of length shows a more significant differ-ence, with the pipe able to encapsulate much longer droplets than thechannel.• The pipe geometry uses the hoop stress to reduce normal stress effects.• A different way in which droplet size may be increased is via introduc-tion of a density difference between fluids. When pumping in the direc-tion of gravity, lighter droplets increase encapsulation volume (heavierif pumping against gravity).• The increase in encapsulation volume due to favourable density dif-ference is not unbounded, but appears to reach a maximum approxi-mately when the frictional pressure gradient of the base flow matchesthe buoyancy, i.e.∂p0∂x= −Byy,0= χ ⇒τˆYgˆRˆyy,0= ρˆ− ρˆd. (7.1)On decreasing the buoyancy beyond this limit, buoyancy stresses con-tribute to breaking of the plug. Simplistically, the main mechanismof buoyancy at the limit (7.1) is to compensate within the droplet forpressure variations in the yielded layer. This then reduces the need fornormal stress gradients between the droplet interface and yield surface.• The computational solutions reveal rather complex stress distributionsinside the plug, not easily inferred from the asymptotic analysis.• Although we have computed elliptical shaped droplets in this paper,in principle a yield stress allows any shape to be encapsulated. Figure7.1 shows some exotic droplet shapes that do not yield the plug.617.2. Limitations• The main results of this thesis have been submitted to Journal of FluidMechanics for publication: “Macro-size drop encapsulation, A. Maleki,S. Hormozi, A. Roustaie, I.A. Frigaard”.7.2 Limitations• The asymptotic solution predicts only the velocity and stress distri-butions outside the plug, but not the stress distribution inside theplug. Although we attempted to derive predictions of the (indetermi-nate) stress fields inside the plug, in order to predict breaking, we wereunsuccessful. Details of unsuccessful attempts are listed in appendix§B.• The asymptotic solution is unable to predict other flow features suchas spacing between successive droplets.• As we assume a priori that the droplet translates rigidly in the plug,there is no way to include the density difference in the asymptotic solu-tion, i.e. we solve a 1 fluid problem with a perturbed droplet boundary.• We have not considered surface tension effects in this thesis. De-pending on the chosen droplet fluid there may anyway be no surfacetension. We note that in using miscible fluids in many practical situa-tions Pe´clet numbers (based on molecular diffusion) are very large andwith no relative droplet motion dispersion is nullified. If the dropletis immiscible the validity of this neglect rests on the relevant stresses,and we expect our results to retain validity provided thatτˆY σˆTκˆ, (7.2)where σˆT and κˆ are the surface tension coefficient and (a representa-tive) radius of curvature of the droplet. Surface tension values betweenyield stress fluids and other fluids are not well known and reliable mea-surement techniques are still being developed, e.g. [9]. Nevertheless,for values typical of water and light oils, with yield stresses of ∼ 10 Pacapillary effects are insignificant on lengthscales & 2 mm, which rep-resents many industrial scale pumping operations. This is not to saythat interesting flow effects might not be produced by encapsulatingimmiscible fluid droplets. Indeed, these may have some advantagesin the forming part of such a flow where it is likely that the plug isunyielded as a droplet is introduced to the main flow.627.3. Future directions7.3 Future directionsA number of areas could be studied within this general direction.• A challenging next step in order to validate/verify our results wouldbe to set up macro-encapsulation experiments.• One can consider more realistic models for the yield stress fluids,e.g. Herschel-Bulkley. However, we do not believe that this changes thecharacteristics of encapsulation dramatically from those found here.• Considering other rheological effects such as visco-elasticity as wellas including surface tension which may aid in the establishment ofencapsulation.• For practical purposes it is necessary to study the formation of droplets,either numerically or experimentally. The question here might be howto control droplet shape and frequency in the formation stage.63Bibliography[1] Adachi, K. & Yoshioka, N. 1973 On creeping flow of a visco-plasticfluid past a circular cylinder. Chem. Eng. Sci. 28 (1), 215 – 226.[2] Al-Mohammad, A.M., Alkhaldi, M.H., Al-Mutairi, S.H. & Al-Zahrani, A.A. 2012 Acidizing-induced damage in sandstone injectorwells: laboratory testing and a case history. SPE J. 17 (3), 885–902.[3] Balmforth, N. & Craster, R. 1999 A consistent thin-layer theoryfor Bingham plastics. J. Non-Newt. Fluid Mech. 84, 65–81.[4] Balmforth, N., Frigaard, I.A. & Ovarlez, G. 2014 Yielding tostress: Recent developments in viscoplastic fluid mechanics. Annu. Rev.Fluid Mech. 46, 121–146.[5] Beaulne, M. & Mitsoulis, E. 1997 Creeping motion of a sphere intubes filled with Herschel-Bulkley fluids. J. Non-Newt. Fluid Mech. 72,55–71.[6] Beris, A.N., Tsamopoulos, J.A. & Armstrong, R.C. 1985 Creep-ing motion of a sphere through a Bingham plastic. J. Fluid Mech. 158,219–244.[7] Bhavaraju, S.M., Mashelkar, R.A. & Blanch, H.W. 1978 Bub-ble motion and mass-transfer in non-newtonian fluids .1. single bubblein power law and Bingham fluids. AIChE J. 24 (6), 1063–1070.[8] Blackery, J. & Mitsoulis, E. 1997 Creeping motion of a sphere intubes filled with a Bingham plastic material. J. Non-Newt. Fluid Mech.70, 59–77.[9] Boujlel, J. & Coussot, P. 2013 Measuring the surface tension ofyield stress fluids. Soft Matt. 9 (25), 5898–5908.[10] Chateau, X., Ovarlez, G. & Trung, K.L. 2008 Homogenizationapproach to the behavior of suspensions of noncolloidal particles in yieldstress fluids. J. Rheol 52 (2), 489–506.64Bibliography[11] Cohen, C.E., Ding, D., Quintard, M. & Bazin, B. 2008 Frompore scale to wellbore scale: Impact of geometry on wormhole growthin carbonate acidization. Chem. Eng. Sci. 63 (12), 3088–3099.[12] Cohen, I., Li, H., Hougland, J.L., Mrksich, M. & Nagel, S.R.2001 Using selective withdrawal to coat microparticles. Science 292,265–267.[13] Cox, B.G. 1964 An experimental investigation of the streamlines inviscous fluid expelled from a tube. J. Fluid Mech. 20 (2), 193–200.[14] De Besses, B.D., Magnin, A. & Jay, P. 2003 Viscoplastic flowaround a cylinder in an infinite medium. J. Non-Newt. Fluid Mech.115 (1), 27–49.[15] Dimakopoulos, Y., Pavlidis, M. & Tsamopoulos, J. 2013 Steadybubble rise in Herschel-Bulkley fluids and comparison of predictions viathe augmented Lagrangian method with those via the Papanastasioumodel. J. Non-Newt. Fluid Mech. 200, 34–51.[16] Dubash, N. & Frigaard, I.A. 2004 Conditions for static bubbles inviscoplastic fluids. Phys. Fluids 16, 4319–4330.[17] Dubash, N. & Frigaard, I.A. 2005 Propagation and stopping of airbubbles in carbopol solutions. J. Non-Newt. Fluid Mech. 142, 123–134.[18] Dunbrack, G. 2013 Interfacial effects in visco-plastic lubrication flows.Master’s thesis, University of Britsh Columbia.[19] Duvaut, G. & Lions, J.L. 1976 Inequalities in Mechanics andPhysics. Springer-Verlag.[20] Fortin, M. & Glowinski, R. 1983 Augmented Lagrangian methods:application to the numerical solution of boundary-value problems. Else-vier Science Publishers B.V.[21] Frigaard, I.A. 2001 Super-stable parallel flows of multiple visco-plastic fluids. J. Non-Newt. Fluid Mech. 100 (1-3), 49–76.[22] Frigaard, I.A. & Nouar, C. 2005 On the usage of viscosity regular-isation methods for visco-plastic fluid flow computation. J. Non-Newt.Fluid Mech. 127 (1), 1–26.[23] Frigaard, I.A. & Ryan, D.P. 2004 Flow of a visco-plastic fluid in achannel of slowly varying width. J. Non-Newt. Fluid Mech. 123, 67–83.65Bibliography[24] Furui, K., Burton, R.C., Burkhead, D.W., Abdelmalek, N.A.,Hill, A.D., Zhu, D. & Nozaki, M. 2012a A comprehensive modelof high-rate matrix-acid stimulation for long horizontal wells in car-bonate reservoirs: part I-scaling up core-level acid wormholing to fieldtreatments. SPE J. 17 (1), 271–279.[25] Furui, K., Burton, R.C., Burkhead, D.W., Abdelmalek, N.A.,Hill, A.D., Zhu, D. & Nozaki, M. 2012b A comprehensive model ofhigh-rate matrix-acid stimulation for long horizontal wells in carbonatereservoirs: part II-wellbore/reservoir coupled-flow modeling and fieldapplication. SPE J. 17 (1), 280–291.[26] Ganan-Calvo, A.M. 1998 Generation of steady liquid microthreadsand micron-sized monodisperse sprays in gas streams. Phys. Rev. Lett.80 (2), 285–288.[27] Glowinski, R. 1984 Numerical methods for nonlinear variational prob-lems. Springer-Verlag New York Inc.[28] Glowinski, R. & Wachs, A. 2011 On the numerical simulation ofviscoplastic fluid flow. Handb. Numer. Anal. 16, 483–717.[29] Gref, R., Minamitake, Y., Peracchia, M.T., Trubetskoy, V.,Torchilin, V. & Langer, R. 1994 Biodegradable long-circulatingpolymeric nanospheres. Science 263, 1600–1603.[30] Hillery, A.M., Lloyd, A.W. & Swarbrick, J. 2001 Drug deliveryand targeting for pharmacists and pharmaceutical scientists. Taylor andFrancis Inc.[31] Holenberg, Y., Lavrenteva, O.M., Liberzon, A., Shavit, U. &Nir, A. 2013 PTV and PIV study of the motion of viscous drops inyield stress material. J. Non-Newt. Fluid Mech. 193, 129–143.[32] Hormozi, S. 2011 Multi-layer flows with yield stress fluids. PhD thesis,University of Britsh Columbia.[33] Hormozi, S., Dunbrack, G. & Frigaard, I.A. 2014 Transient be-haviour and shaped interface formation in non-equilibrium multilayerflows. Phys. Fluids p. under publication.[34] Hormozi, S. & Frigaard, I.A. 2012 Nonlinear stability of a visco-plastically lubricated viscoelastic fluid flow. J. Non-Newt. Fluid Mech.169, 61–73.66Bibliography[35] Hormozi, S., Martinez, D. M. & Frigaard, I.A. 2011c Stablecore-annular flows of viscoelastic fluids using the visco-plastic lubrica-tion technique. J. Non-Newt. Fluid Mech. 166 (23-24), 1356–1368.[36] Hormozi, S., Wielage-Burchard, K. & Frigaard, I.A. 2011aMulti-layer channel flows with yield stress fluids. J. Non-Newt. FluidMech. 166 (5-6), 262–278.[37] Hormozi, S., Wielage-Burchard, K. & Frigaard, I.A. 2011bEntry, start up and stability effects in visco-plastically lubricated pipeflows. J. Fluid Mech. 673, 432–467.[38] Huen, C.K., Frigaard, I.A. & Martinez, D.M. 2007 Experimentalstudies of multi-layer flows using a visco-plastic lubricant. J. Non-Newt.Fluid Mech. 142 (1-3), 150–161.[39] Jaworek, A. 2008 Electrostatic micro and nanoencapsulation andelectroemulsification: A brief review. J. Microencapsul. 25 (7), 443–468.[40] Kalfayan, L. 2008 Production enhancement with acid stimulation,2nd edn. Penwell Corportation.[41] Lavrenteva, O.M., Holenberg, Y. & Nir, A. 2009 Motion ofviscous drops in tubes filled with yield stress fluid. Chem. Eng. Sci.64 (22), 4772–4786.[42] Lipscomb, G.G. & Denn, M. 1984 Flow of Bingham fluids in complexgeometries. J. Non-Newt. Fluid Mech. 14, 337–346.[43] Lister, J.R. 1989 Selective withdrawal from a viscous 2-layer system.J. Fluid Mech. 198, 231–254.[44] Liu, B.T., Muller, S.J. & Denn, M.M. 2002 Convergence of aregularization method for creeping flow of a Bingham material about arigid sphere. J. Non-Newt. Fluid Mech. 102, 179–191.[45] Liu, B. T., Muller, S. J. & Denn, M. M. 2003 Interactions oftwo rigid spheres translating collinearly in creeping flow in a Binghammaterial. J. Non-Newt. Fluid Mech. 113 (1), 49–67.[46] Liu, X.H., Ormond, A., Bartko, K., Li, Y. & Ortoleva, P.1997 A geochemical reaction-transport simulator for matrix acidizinganalysis and design. J. Petrol. Sci. Eng. 17 (1-2), 181–196.67Bibliography[47] Loscertales, I.G., Barrero, A., Guerrero, I., Cortijo, R.,Marquez, M. & Ganan-Calvo, A.M. 2002 Micro/nano encapsuta-tion via electrified coaxial liquid jets. Science 295, 1695–1698.[48] Mahaut, F., Chateau, X., Coussot, P. & Ovarlez, G. 2008Yield stress and elastic modulus of suspensions of noncolloidal particlesin yield stress fluids. J. Rheol. 52 (1), 287–313.[49] Maheshwari, P. & Balakotaiah, V. 2013 Comparison of carbon-ate HCl acidizing experiments with 3D simulations. SPE Prod. Oper.28 (4), 402–413.[50] Merkak, O., Jossic, L. & Magnin, A. 2009 Migration and sedimen-tation of spherical particles in a yield stress fluid flowing in a horizontalcylindrical pipe. AIChE J. 55 (10), 2515–2525.[51] Moyers-Gonzalez, M.A., Frigaard, I.A. & Nouar, C. 2004 Non-linear stability of a visco-plastically lubricated viscous shear flow. J.Fluid Mech. 506, 117–146.[52] Moyers-Gonzalez, M.A., Frigaard, I.A. & Nouar, C. 2010 Sta-ble two-layer flows at all Re; visco-plastic lubrication of shear-thinningand viscoelastic fluids. J. Non-Newt. Fluid Mech. 165 (23-24), 1578–1587.[53] O’Donovan, E.J. & Tanner, R.I. 1984 Numerical study of the Bing-ham squeeze film problem. J. Non-Newt. Fluid Mech. 15 (1), 75–83.[54] Ovarlez, G., Bertrand, F., Coussot, P. & Chateau, X. 2012Shear-induced sedimentation in yield stress fluids. J. Non-Newt. FluidMech. 177, 19–28.[55] Ovarlez, G., Bertrand, F. & Rodts, S. 2006 Local determinationof the constitutive law of a dense suspension of noncolloidal particlesthrough magnetic resonance imaging. J. Rheol. 50 (3), 259–292.[56] Papanastasiou, T.C. 1987 Flows of materials with yield. J. Rheol.31 (5), 385–404.[57] Potapov, A., Spivak, R., Lavrenteva, O.M. & Nir, A. 2006Motion and deformation of drops in Bingham fluid. Ind. Eng. Chem.Res. 45 (21), 6985–6995.68Bibliography[58] Priya, A. R. Sathiya, Muralidharan, V. S. & Subramania, A.2008 Development of novel acidizing inhibitors for carbon steel corro-sion in 15% boiling hydrochloric acid. Corros. Sci. 64 (6), 541–552.[59] Putz, A., Burghelea, T.I., Frigaard, I.A. & Martinez, D.M.2008 Settling of an isolated spherical particle in a yield stress shearthinning fluid. Phys. Fluids 20 (3).[60] Putz, A. & Frigaard, I.A. 2010 Creeping flow around particles ina Bingham fluid. J. Non-Newt. Fluid Mech. 165, 263–280.[61] Putz, A., Frigaard, I.A. & Martinez, D.M. 2009 On the lubri-cation paradox and the use of regularisation methods for lubricationflows. J. Non-Newt. Fluid Mech. pp. 62–77.[62] Quraishi, M.A. & Jamal, D. 2000 Technical note: CAHMT - A newand eco-friendly acidizing corrosion inhibitor. Corros. Sci. 56 (10), 983–985.[63] Ratnakar, R.R., Kalia N. & Balakotaiah, V. 2013 Modeling,analysis and simulation of wormhole formation in carbonate rocks within situ cross-linked acids. Chem. Eng. Sci. 90 (0), 179–199.[64] Roquet, N. & Saramito, P. 2003 An adaptive finite element methodfor bingham fluid flows around a cylinder. Comput. Methods in Appl.Mech. Eng. 192, 3317–3341.[65] Roustaei, A. & Frigaard, I.A. 2013 The occurrence of fouling layersin the flow of a yield stress fluid along a wavy-walled channel. J. Non-Newt. Fluid Mech. 198, 109–124.[66] Ryan, D.P. 2004 Slow flow of bingham fluid in a gap of slowly varyingwidth. Master’s thesis, University of Britsh Columbia.[67] Saramito, P. & Roquet, N. 2001 An adaptive finite element methodfor viscoplastic fluid flows in pipes. Comput. Methods in Appl. Mech.Eng. 190 (40-41), 5391–5412.[68] Singh, J.P. & Denn, M.M. 2008 Interacting two-dimensional bubblesand droplets in a yield-stress fluid. Phys. Fluids 20 (4).[69] Smith, C.F. & Hendrick, A.R. 1965 Hyrofluoric acid stimulation ofsandstone reservoirs. J. Petrol. Technol. 17 (2), 215–222.69Bibliography[70] Stein, S & Buggisch, H 2000 Rise of pulsating bubbles in fluids witha yield stress. J. Appl. Math. Mech. 80 (11-12), 827–834.[71] Tokpavi, D.L., Magnin, A. & Jay, P. 2008 Very slow flow of Bing-ham viscoplastic fluid around a circular cylinder. J. Non-Newt. FluidMech. 154 (1), 65–76.[72] Walton, I.C. & Bittleston, S.H. 1991 The axial-flow of a Binghamplastic in a narrow eccentric annulus. J. Non-Newt. Fluid Mech. 222,39–60.[73] Windbergs, M., Zhao, Y., Heyman, J. & Weitz, D.A. 2013Biodegradable core-shell carriers for simultaneous encapsulation of syn-ergistic actives. J. Am. Chem. Soc. 135 (21), 7933–7937.[74] Zhao, Y., Shum, H.C., Adams, L.A., Sun, B., Holtze, C., Gu, Z.& Weitz, D.A. 2011 Enhanced encapsulation of actives in self-sealingmicrocapsules by precipitation in capsule shells. Langmuir 27 (23),13988–13991.[75] Zou, C.J., Liao, W.J., Zhang, L. & Chen, H.M. 2011 Study onacidizing effect of beta-cyclodextrin-PBTCA inclusion compound withsandstone. J. Petrol. Sci. Eng. 77 (2), 219–225.[76] Zuidam, N.J. & Nedovi, V.A. 2010 Encapsulation technologies foractive food ingredients and food processing . Springer.70Appendix AAlgebraic DetailsA.1 Channel Poiseuille flow of a Bingham fluidHere we want to study the plane Poiseuille solution. We assume the flowis fully developed and there is no x-gradient. We also neglect the gravityhere assuming that the channel is horizontal. Due to the symmetry of theproblem, without of loss of generality, we can only consider y > 0. Finally,in order to find pressure gradient we take flow rate to be unity. With theseassumptions 2.9b simplifies to−∂p∂x+∂τxy∂y= 0 (A.1)Knowing that p = p(x) and the stress has to be zero at the centerline dueto symmetry, we write:τxy = −∂p∂xy (A.2)Yielding happens where τ = |τxy| ≥ B and therefore position of plug in-terface is characterised by B = −yy,0∂p∂x . In the yielded region y ≥ yy,0,∂UP∂y < 0 and therefore constitutive relations simplify to τxy =∂UP∂y − B.Thus∂UP∂y=∂p∂xy −B =∂p∂x(y − yy,0) ⇒ UP =∂p∂x((y − yy,0)2 − (1− yy,0)2)(A.3)where we satisfy no slip condition at the wall. This can be rearranged like:UP (y) =up,0 if |y| ≤ yy,0,up,0[1−(|y| − yy,0)2(1− yy,0)2]if |y| > yy,0,(A.4)in which up,0 is the plug velocity and is given byup,0 =B2yy,0(1− yy,0)2, (A.5)71A.1. Channel Poiseuille flow of a Bingham fluidThe position of the yield surface is found by demanding a flow rate ofunity:∫ 10UPdy =∫ yy,00up,0dy +∫ 1yy,0UP (y)dy = 1 (A.6)which leads to Buckingham equation:f(y) := y3 − (3 +6B)y + 2 = 0, and f(yy,0) = 0 (A.7)One can notice that for any B > 0, f(0).f(1) < 0 which means∃yy ∈ (0, 1) : f(yy) = 0 (A.8)regardless of the value of B. Also f ′(y) < 0 for y ∈ (0, 1), so that there maybe only a single root.In limit when B −→ 0, we can write yy = yy0 +Byy1 +B2yy2. Pluggingthis expansion in the Buckingham Equation, collecting powers of B andequating them to zero yields:6Byy0 = 0 ⇒ yy0 = 0 (A.9a)y3y0 − 3yy0 − 6yy1 + 2 = 0 ⇒ yy1 =13(A.9b)B(3y2y0yy1 − 3yy1 − 6yy2)= 0 ⇒ yy2 = −16(A.9c)and finally we can write:B −→ 0 yy ∼B3−B26(A.10)This is consistent with our expectations. In particular, for a Newtonian fluid(B = 0), the thickness of plug region is zero (yy = 0).On the other side of the spectrum we consider the case of B −→∞. Thezero order solution corresponds to the blockage of channel when the channelin entirely filled with plug region yy = 1. Looking for the perturbationsaround this root, we can expand yy = 1 + αBp; where p and α have to befound. Plugging into the Buckingham Equation:3αBp + 3α2B2p −6B− 3αBp − 6αBp−1 +O(3p) = 0 (A.11)Given that α 6= 0, the only possible value for p is −1/2. Then we can findα = −√2 and therefore,B −→∞ yy ∼ 1−√2B(A.12)72A.2. Asymptotic form of Bingham constitutive relationsA.2 Asymptotic form of Bingham constitutiverelationsA.2.1 channel geometryIn this section we derive the asymptotic form of Bingham constitutive rela-tions in Cartesian geometry. We adopt a regular asymptotic expansion inthe form of(u, v, p) = (u0, v0, p0) + δ(u1, v1, p0) + δ2(u2, v2, p2) + · · · (A.13)where u and v represent dimensionless components of velocity and p is di-mensionless pressure introduced in §3.1. 7 We haveγ˙xx = 2∂u∂x= 2∂u0∂x+ δ2∂u1∂x+ δ22∂u2∂x+O(δ3) (A.14a)γ˙xy =∂u∂y+ δ2∂v∂x=∂u0∂y+ δ∂u1∂y+ δ2(∂u2∂y+∂v0∂x)+O(δ3) (A.14b)Substituting (A.14) in (1.4c) yields:γ˙ =(∂u0∂y2+ δ2∂u0∂y∂u1∂y+ δ2[∂u1∂y2+ 2∂u0∂y(∂u2∂y+∂v0∂x)+ 4∂u0∂x2])12=∣∣∣∣∂u0∂y∣∣∣∣1 + δ2∂u1∂y∂u0∂y+ δ2∂u1∂y2+ 2∂u0∂y(∂u2∂y +∂v0∂x)+ 4∂u0∂x2∂u0∂y212(A.15)Knowing that δ  1, we can use Taylor expansion to find:γ˙−1 =1∣∣∣∂u0∂y∣∣∣1− δ∂u1∂y∂u0∂y− δ24∂u1∂y2+ 8∂u0∂y(∂u2∂y + 4∂v0∂x)+ 16∂u0∂x2+ 3∂u1∂y28∂u0∂y2(A.16)Now we haveτxy =(1 +Bγ˙)γ˙xy =∂u0∂y+B sgn(∂u0∂y) + δ∂u1∂y+O(δ2) (A.17a)τxx =(1 +Bγ˙)γ˙xx = 21 +B∣∣∣∂u0∂y∣∣∣∂u0∂x+O(δ) (A.17b)This analysis is taken from [66] with minor changes.7The uppercase letters of §3.1 are shown in lower case here.73A.2. Asymptotic form of Bingham constitutive relationsA.2.2 pipe geometryIn this section we derive the asymptotic form of Bingham constitutive rela-tions in cylindrical geometry. We adopt a regular asymptotic expansion tothe form of(u,w, p) = (u0, w0, p0) + δ(u1, w1, p0) + δ2(u2, w2, p2) + · · · (A.18)where u and w represent dimensionless radial and axial velocities respec-tively and p is dimensionless pressure all introduced in §6.2.The uppercaseletters of §6.2 are shown in lower case here. Using the identity:γ˙rr + γ˙θθ + γ˙zz = 0,we rewrite (6.10b) as:γ˙ =√γ˙2rz − δ2(γ˙rrγ˙θθ + γ˙rrγ˙zz + γ˙rrγ˙zz)(A.19)In addition, we knowγ˙rz =∂u0∂r+ δ∂u1∂r+ δ2(∂u2∂r+∂v0∂z) +O(δ3) (A.20a)γ˙rr = 2∂v0∂r+O(δ) (A.20b)γ˙θθ = 2v0r+O(δ) (A.20c)γ˙zz = 2∂u0∂z+O(δ) (A.20d)To go on, we need to calculate 1γ˙ , so we haveγ˙ =∣∣∣∣∂u0∂r∣∣∣∣[1 + δ(2∂u1∂r∂u0∂r)+ δ2(2∂u0∂r (∂u2∂r +∂v0∂z ) +∂u1∂r2− 4(v0r∂v0∂r +∂v0∂r∂u0∂z +v0r∂u0∂z )∂u0∂r2)]1/2(A.21)which results inγ˙−1 =1∣∣∣∂u0∂r∣∣∣[1− δ∂u1∂r∂u0∂r− δ2(∂u0∂r (∂u2∂z +∂v0∂z )− 2(v0z∂v0∂r +∂u2∂z∂v0∂r +v0r∂u0∂z )−52∂u1∂r2∂u0∂r2)](A.22)74A.3. Perturbations are O(δ2)Thusτrz =∂u0∂r+B sgn(∂u0∂r) + δ∂u1∂r+O(δ2) (A.23a)τrr =2∂v0∂r(1 +B|∂u0∂r |)+O(δ) (A.23b)τθθ =2v0r(1 +B|∂u0∂r |)+O(δ) (A.23c)τzz =2∂u0∂z(1 +B|∂u0∂r |)+O(δ) (A.23d)A.3 Perturbations are O(δ2)A.3.1 yT − yy,0 ∼ O(δ2)we define ξ1 the departure of 0-order plug interface (yy) from its Poiseuillecounterpart : yy = yy,0 + δξ1. Substituting this into (3.13) gives:O(1) : y3y,0 − yy,0(3 + 6/B) + 2 = 0, (A.24a)O(δ) : 3y2y,0ξ1 − ξ1(3 + 6/B) + 6h¯B= 0, =⇒ ξ1 =h¯1 + 0.5(1− y2y,0)(A.24b)The first expression is simply the Buckingham equation (2.16), as expected.The second expression defines ξ1. Now, since up is defined by yy and h, westraightforwardly expand 3.12 in powers of δ and writeη =Bδ(up,0 − up(X)) =B2δ((1− yy,0)2yy,0−(1− (yy,0 + δξ1))2(yy,0 + δξ1)− δh¯)=B2δ((1− yy,0)2yy,0−(1− yy,0)2yy,0− δ(1− yy,0)(ξ1 − h¯− 3yy,0ξ1 + yy,0h¯y2y,0)=B2y2y,0((1− yy,0)(ξ1 − h¯+ yy,0ξ1 + yy,0h¯))+O(δ)(A.25)75A.3. Perturbations are O(δ2)This confirms that η = O(1) as has been claimed earlier. To evaluate yT ,we need to find ∂U1∂y (X, yy) and(∂2U0∂Y 2 (X, yy))−1. From (3.19) we have:∂U1∂y(X, yy) = −2ηyy + 2(yy − 1)2= −2η(yy,0 + δξ1) + 2((yy,0 + δξ1)− 1)2= −2η(yy,0 + 2(yy,0 − 1)2− δξ1yy,0 + 5(ξ1 − 1)3)+O(δ2)(A.26)and from (3.11)(∂2U0∂Y 2(X, yy))−1= −yy − δh¯B= −yy,0B+O(δ) (A.27)Substituting the expansion of yy into (3.20) and incorporating (A.25), (A.26)and (A.27), we have:yT−yy,0 = δ(ξ1−∂U1∂y (X, yy)∂2U0∂Y 2 (X, yy)= δ(ξ1 −yy,0 + 2yy,0 − 1yy,0ξ1 + yy,0h¯+ yy1 − h¯yy,0)+O(δ2)(A.28)and then with A.24b:yT − yy,0 =δh¯By4y,0 +By3y,0 − 3By2y,0 − 6y2y,0 −Byy,0 − 6yy,0 + 2Byy,0(yy,0 − 1)(By2y,0 −B − 2)+O(δ2)=hB(yy,0 + 1)(y3y,0 − 3yy,0 −6Byy,0 + 2)yy,0(yy,0 − 1)(By2y,0 −B − 2)+O(δ2)(A.29)Finally, taking Buckingham equation (2.16) into account, we seeyT − yy,0 ∼ O(δ2) (A.30)A.3.2 U0 + δU1 = UP +O(δ2)Similar to the previous section, we define ξ2 the departure of 0-order velocityfrom poiseuille flow velocity: U0 = UP + δξ2. Our goal is to show thatξ2 + U1 ∼ O(δ).Starting from 3.19 and expanding yy = yy,0 + δξ1 we have:U1(X,Y ) = −η(Y − 1)(3Y + 3Y yy,0 − 1− 4y2y,0 − yy,0)(yy,0 − 1)3+O(δ) (A.31)76A.3. Perturbations are O(δ2)and then substituting η from (A.25) :U1(X,Y ) = Bh¯(Y − 1)(3Y + 3Y yy,0 − 1− 4y2y,0 − yy,0)(By2y,0 − 2yy,0B − 2yy,0 +B)(yy,0 − 1)2y2y,0(−2−B +By2y,0)(A.32)From definition of ξ2, we can findξ2 =1δ(U0 − UP ) =B2y2y,0(Y − 1)(ξ1Y − Y h¯+ 2h¯yy,0 + ξ1 − h¯) (A.33)After substitution of ξ1 from 3.21b and algebraic simplification we can showξ2 + U1 = α(y3y,0 − yy,0(3 +6B) + 2)+O(δ) (A.34)in which α is a collection of numerous terms multiplied to each other. Thisthen leads to the conclusion:ξ2 + U1 ∼ O(δ) =⇒ U0 + δU1 = UP +O(δ2) (A.35)77Appendix BUnsuccessful attemptsWe spent a lot of effort in order to find a mathematical relation whichestimates the maximum size of encapsulated drop. Although the results werenot successful, it is still worth mentioning the approach and the challengesthat we faced, hoping that could be a guideline for further studies.B.1 Method of islands of broken plug of [23]In [23], the authors gave an estimate to the maximum of amplitude of awavy channel which does not break the plug region. They performed thisanalysis by balancing continuity and momentum equations on the islandsof the broken plug. We took the same approach here and find the averagevalue of extensional stress along any vertical line inside the channel andsupposedly, this has to be B/δ so that the plug yields.With rearranging and combing the continuity and x-momentum equa-tions together, and then using Green’s theory one can straightforwardlyshow ∮∂AJ1dy + J2dx = 0 (B.1)whereJ1 =− p+ δ2τxx − δReu2 (B.2a)J2 =− τxy + δReuv (B.2b)and A is any simple closed piecewise smooth curved in x − y plane. Wechoose A as shown in Figure B.1 The integration in B.1 can be divided intofour integration on four segments of the boundary (as shown in Figure ).We know yT = yy,0 +O(δ2), we can replace yT by yy,0 to make the analysissimpler. Therefore we have∮(1)J1dy +∮(2)J2dx+∮(3)J1dy +∮(4)J1dy + J2dx = 0 (B.3)78B.1. Method of islands of broken plug of [23]Figure B.1: Methods of island of broken plug∮(1)J1dy =∫ yy,0h(xb)J1(xb, y)dy (B.4a)∮(2)J2dx =∫ 0xbJ2(x, yy,0)dx (B.4b)∮(3)J1dy =∫ yy,0h(0)J1(0, y)dy (B.4c)∮(4)J1dy + J2dx =∫ xb0(h′J1(x, h(x))+ J2(x, h(x)))dx (B.4d)On y = yy,0, we know that τxy = B and v = 0, so∫ 0xbJ2(x, yy,0)dx =∫ 0xb−Bdx = Bxb (B.5)At x = 0, due to symmetry, τxx(0, y) = 0 and p = 0 which gives∫ h(0)yy,0J1(0, y)dy = −δReu2p(h(0)− yy,0)(B.6)On the boundary of drop (y = h(x)), we know τxy = 0 and v = 0, so∫ xb0J2(x, h(x))dx = 0 (B.7)79B.2. Finding an admissible form of stressand using integration by parts∫ xb0h′J1 (x, h(x)) dx = hJ1 (x, h(x))]xb0−∫ xb0h(x)∂J1∂xdx = −δReu2p (h(xb)− h (0))(B.8)Inside the plug region we know horizontal velocity in plug region is constantu(xb, y) = up and vertical velocity, up to O(δ), is zero. Therefore,∫ yy,0h(xb)J1(xb, y)dy = −δReu2p (yy,0 − h(xb))−∫ yy,0h(xb)pdy+δ2∫ yy,0h(xb)τxx(xb, y)dy(B.9)substitution givesδ2∫ yy,0h(xb)τxx(xb, y)dy =∫ yy,0h(xb)pdy −Bxb (B.10)We define I := δ∫ yy,0h(xb)τxx(xb, y)dy and compare it with B as a measureto determine whether drop breaks the plug or not. To proceed we need toguess the variation of the normal stresses inside the channel. The pressurechanges from Ppois = − Byy,0x on the surface of plug region to 0 on the dropinterface, assuming that there is no density difference. As an example, if weadopt a linear profile for pressure, and look for geometrical conditions whichsatisfy I ≤ B, we find 3L2δ ≤ 1. This result is not intuitively consistent withour earlier results.B.2 Finding an admissible form of stressInside the plug region the Bingham constitutive relations do not fully deter-mine the stress behaviour. This means that any form of stress distribution,provided that it satisfies Stokes momentum equations, is admissible.One idea that could lead to find an estimation of maximum size of dropsis to find an admissible solution and play with the geometry of the dropuntil the total stress exceeds the yield stress (B).We start with choosing a linear profile for τxy inside the plug regionwhich satisfies the two boundary conditions that we have at the interfaceof drop and plug region(τxy = −By−h(x)yy0−h(x)). Using horizontal and verticalmomentum equation, we can write∂τxy∂y= −B1yy0 − h(x)= −∂σxx∂x=∂p∂x−∂τxx∂x(B.11)80B.2. Finding an admissible form of stressand∂τxy∂x= −B(yy0 − y)(yy0 − h(x))2h′(x) =∂p∂y−∂τyy∂y=∂p∂y+∂τxx∂y(B.12)and then2∂2τxx∂x∂y= B∂∂x[ yy0 − y(yy0 − h(x))2h′(x)](B.13)Knowing that at x = 0 τxx = 0 and h′(x) = 0τxx = −B4( yy0 − yyy0 − h(x))2h′(x) + f(x) (B.14)From B.12 we can conclude thatp(x, y) = −B4( yy0 − yyy0 − h(x))2h′(x) + g(x) (B.15)which in turn results inσyy = −p+ τyy =B2( yy0 − yyy0 − h(x))2h′(x)− f(x)− g(x) (B.16)we know that at y = h(x), σyy = 0 and therefore f(x) + g(x) = B2 h′(x).On the other hand, at y = yy0, σyy = −p0(x) = − Byy0x and thereforef(x) + g(x) = Byy0x. It is not possible to satisfy both these two boundaryconditions simultaneously.An alternative method is to propose a different function for τxy in theform of τxy = −Bf(η) where η =y−h(x)yy0−h(x)with the conditions of f(0) = 0and f(1) = 1; e.g. τxy = −12B(η2 + η). However, this method is alsoproblematic as the integration required to find τxx makes the equationsalmost unsolvable analytically.81


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



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"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items