Ind. Eng. Chem. Res. 2006, 45, 6985-6995
6985
Motion and Deformation of Drops in Bingham Fluid Alexander Potapov, Roman Spivak, Olga M. Lavrenteva, and Avinoam Nir* Department of Chemical Engineering, TechnionsIsrael Institute of Technology, Haifa 32000, Israel
The axisymmetric dynamics of single and interacting deformable fluid inclusions in yield stress fluid are investigated. The drops move in a cylindrical tube under the action of gravity. The velocities of the drops, deformation patterns, and the evolution of yield domains were examined for a variety of physical parameters including the interfacial tension, the yield stress, and the geometry of the system. Particular attention is given to the evolution of the yielded regions in the ambient matrix and to the deformation of the interfaces of the inclusions. The results of simulations show that, for a single drop, a quasistationary state is established after a relatively long transient period. The interaction of two drops results in a substantial extension of the yielded domain and in the acceleration of the motion. Two equal-size drops tend to approach each other and coalesce. Afterward, a single large fluid particle is formed that sometimes breaks up forming one large leading drop and one or several small trailing satellites. 1. Introduction Yield stress, or viscoplastic, materials constitute an important class of non-Newtonian materials. They behave like solids when they are exposed to low levels of stress. They might remain undeformed or respond to the imposed load as an elastic or viscoelastic solid and exhibit little deformation. If the stress level exceeds a critical value, commonly known as the yield stress, the material yields and flows as a non-Newtonian fluid. Thus, in any particular setup, there may be regions where part of the system reacts as a solid and other parts flow under the local stress. These regions are known as unyielded and yielded, respectively, and they are separated by a sharp interface the location of which is a priori unknown. Hence, when a pressure difference or a shear stress is applied on a volume of a yield stress fluid in a duct, regions of higher stress level may flow while other regions may move as solid plugs transported by the yielded flowing material, depending on the proximity to the duct walls. Similarly, when a body moves in a yield stress material under gravity or the action of another force field, regions of flowing fluid appear near the body surface while total stagnation is experienced elsewhere. Evidently, such flow patterns are completely dissimilar to those exhibited by any other non-Newtonian material and they have great importance since they govern the ability to process these fluids as well as dominate transport phenomena in them. The most celebrated yield stress material is the Bingham fluid. Known yield stress fluids are encountered in almost all aspects of life. They include mud, cements, lava, glues and paints, various foodstuffs, fermentation broths, foams, suspensions, emulsions, gels, rocket fuels, and other polymer mixtures. Hence, they are of interest because they appear in various natural phenomena as well as a variety of technological applications, these ranging from biotechnology to civil and chemical engineering, from space industry to military applications, and from food and cosmetic production to sophisticated injectable and transdermal drug delivery systems. An assortment of fluids with a yield stress is considered in Bird et al.1 Another interesting review of yield stress materials and associated phenomena, though with the intention of challenging the yield stress characterization, is given by Barnes.2 * To whom correspondence should be addressed. Tel.: +972 829 1921. Fax: +972 829 7256. E-mail:
[email protected].
The rich variety of systems and processes that involve yield stress fluids spans an extremely wide range of scales and applications, from the giant flow of lava or macroscopic motion of mud slurries (Coussot3) down to the nanometric scale of gel coated artificial vesicles in advanced cardiovascular drug delivery practices in micro- and nanoemulsions (see, e.g., a comprehensive collection of reviews by Nielloud and MartiMestres4). In the chemical industry, they are found as suspensions and emulsions, where interesting relations between surface chemistry and rheology are found (see various works by Boger and co-workers, e.g., refs 5 and 6). In biotechnology, they result as fermentation broths in bioreactors. Our toothpaste is a yield stress material. Various cosmetic commodities such as creams and hair gels are members of this family. Medication ointments stored in tubes do not flow unless squeezed. Many processed foods exhibit this behavior as well. Mayonnaise, peanut butter, jams, and ketchup are just a few to be mentioned. Various polymer blends are yield stress materials. Thus, energetic materials as solid rocket fuel behave as such when manufactured and processed (Kang et al.7). On the small scale appear devices that are developed for targeted oral, transdermal, and cardiovascular drag delivery in the form of emulsions that are comprised of compound nanocapsules with a gelatinous biodegradable envelope (e.g., Gallardo et al.8). The analysis of the mechanics of yield stress materials has been of continuous interest since the pioneering work of Bingham,9 and several invariant models for the rheology were formulated throughout the years. Oldroyd10,11 and Prager12 employed a von Mises yield criterion for the onset of flow beyond a critical stress. A popular generalization of the Bingham constitutive description is known as the Herschel-Bulkley model (see, e.g., Burgos et al.13). Very few analytic solutions have been obtained for the flow of yield stress fluids. The reason is, of course, the discontinuous and nonlinear nature of the Bingham constitutive model. Only very simple geometries allow a determination of the location of the interface between the yielded and the unyielded regions as well as satisfaction of all the necessary compatibility conditions there. Such solutions are naturally limited to rectilinear and unidirectional lubrication flows (see Bird et al.1), and some of them contain inaccuracies (Lipscomb and Denn14). Thus, much effort was devoted to the development of variational principles and variational methods for creeping flows (Prager,12
10.1021/ie051222e CCC: $33.50 © 2006 American Chemical Society Published on Web 04/07/2006
6986
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006
Mosolov and Miasnikov,15 Yoshioka and Adachi,16 Duvaut and Lions,17 Glowinski,18 and more recent works by Huilgol, see, e.g., ref 19 and a series of works by Frigaard and co-workers20). These methods are most effective in mathematical analysis of the problem (existence and uniqueness theorems). Also, on the basis of volume and surface integrals, the variational unequalities give upper and lower bounds and good estimates for integral quantities such as drag force. However, they are not very useful to describe pointwise characteristics of the flow such as velocity and stress fields. Thus, ever since the analysis of Glowinski et al.21 that combined inequalities with numerical methods and minimization techniques to calculate viscoplastic flows, the focus on the solutions of such flows was via numerical techniques. Various modifications of finite difference and finite element methods were used extensively, in particular for solving for the dynamics of inclusions in yield stress materials. Some of these works are reviewed below in conjunction with the motion of inclusions in such domains. Dispersed phases are common in yield stress materials. Often the interaction between the inclusions renders the material a yield stress fluid, and it is the micro- or nanoscale composition that is responsible for the special behavior, i.e., providing the material with the ability to behave like a solid when exposed to low stress levels yet flow like a fluid when the forcing is increased. Understanding the dynamics of the minute scale constituents is thus most important to understanding the properties, behavior, and processing ability of these materials. In some cases, the dispersion is a natural result of the process or the manufacturing sequence. Such are red mud suspensions in aluminum processing waste streams (Nguyen and Boger22,23) or phase-separated block copolymer melts (Spaans and Williams24). In other cases, the inclusions are mixed into the system as favorable additives to be thoroughly dispersed. Economic reasons drive the stuffing of various cosmetic gels with tiny insoluble air bubbles since it is a commodity sold by volume. Similarly, the difference of taste between aerated and solid chocolate is a driving force in the food industry. Most often, though, the dispersed entities are an undesired phase that has to be controlled if not separated and removed. Unexpected flow of high pressure gas from rock formations into the drilling mud (a viscoplastic material), and its rise into the well bore, may result in gas kicks that can lead to blowouts at the surface with consequent environmental damage and risk of life (Santos and Azar25). In poured rocket fuels, that are basically suspensions of small inclusions in a continuous polymeric matrix, the particles tend to migrate, and the fuel may lose homogeneity which can result in uneven burning and propulsion. Furthermore, gas bubbles trapped in the fuel can become focal points that develop into hot spots due to rapid compression and can cause self-ignition of the explosive material (Kang et al.7). Such examples, and others, indicate that a basic understanding of the dynamics of small solid and fluid inclusions in yield stress materials is of utmost importance. Most of the studies of the motion of inclusions in viscoplastic materials were devoted so far to spherical solid particles. Early works (e.g., Andres,26 Valentic and Whitmore,27 and Yoshioka et al.28) were on a sphere in an infinite medium, while more recent publications (Beris et al.,32 Blackery and Mitsoulis,29 Beauline and Mitsoulis,30 and Liu et al.31) dealt with numerical simulation of such dynamics as well as the motion inside finite tubes. In an early study,27 a crude approximate theory was developed based on the assumption of a spherical shape of the yielded region. More complicated shapes of this region, offered
by the calculation, were considered by Beris et al.32 In particular, they reported that small unyielded solid regions should be present near the stagnation points at the front and rear poles of the sphere. Beauline and Mitsoulis30 predicted the existence of “unyielded islands” at the equatorial plane of the sphere. However, as argued recently,33 such domains cannot exist at the plane of symmetry for a translation of a single sphere. Such arguments were raised by others in earlier publications.14,31,32 Furthermore, Smyrnaios and Tsamopoulos34 have also recently demonstrated that an unyielded region cannot exist at the plane of symmetry in squeeze flow. Thus, inaccurate stress prediction depends on the particular smoothed rheological model used in the numerical scheme. Several experimental investigations of particles falling in viscoplastic media (Ansley and Smith,35 de Plessis and Ansley,36 Ito and Kajachi,37 and Atapattu et al.38) were focused on measuring the drag coefficient as a function of the liquid viscosity and yield stress. Recently, these studies were extended to nonspherical solid inclusions (Jossic and Magnin39). More recent approaches use methods that follow the augmented Lagrangian approach developed originally by Glowinski,18,21 and they do seem to resolve the shape of the unyielded regions correctly (Roquet and Saramito33). The interaction between two rigid spheres translating with equal velocities in a Bingham material was studied recently by Liu et al.40 They report two striking features that emerge from the calculations: the particles exhibit negligible interaction unless they are in close proximity (separation of less than four sphere radii); when in close proximity, there are some complex regions of unyielded and yielded domain separating the two spheres. When the spheres are approaching each other, the unyielded region separating them diminishes except perhaps the cap zone reported for isolated particles.14,15 Recent experimental studies of a laminar flow of a suspension of spheres in a gelled fluid through a sudden expansion (Jossic and Magnin41,42) revealed migration and segregation of inclusions that can be employed to create organized structures. No theoretical studies of this interesting phenomenon are available so far. The literature on the motion of drops and bubbles in viscoplastic fluids is so far quite limited. One of the early studies of fluid particle dynamics in yield stress materials was reported by Bhavaraju et al.,43 who employed perturbation in the Bingham parameter to study the motion of a nondeformable bubble in an unbounded medium. Stein and Buggisch44 applied perturbation techniques to demonstrate that bubbles trapped in a viscoplastic medium can be made to rise by applying an external oscillatory pressure. Experimental results confirming these predictions are reported there as well. A recent work by Dubash and Frigaard45 addresses the propagation of large bubbles in cylindrical columns and uses variational principles (Prager12) to establish bubble stopping conditions in the tube. Li and Renardy46 reported on numerical simulations of deformation and rupture of a viscous drop in a completely yielded matrix of a viscoplastic fluid sheared between two solid plates. There are no studies reporting the motion of deformable inclusions in Bingham type fluids in the presence of unyielded regions. It should be pointed out that, while it is well-known that small drops and bubbles do not deform when slowly translating in a Newtonian viscous fluid unless other effects such as inertia or the proximity of walls are involved, this is not necessarily the general case for translation in non-Newtonian media except when the interfacial tension is extremely large. Thus, a deformation should be expected. Furthermore, the shape of the yielded domains surrounding drops or bubbles is not necessarily similar
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006 6987
to that found for solid spheres (Liu et al.31). It can be anticipated that the free flow of the fluid-fluid interface would prevent the accumulation of solidified caps near the stagnation points at the poles. In this communication, we report a study of singleand two-drop deformation and interaction when falling freely under gravity in a confined vessel. For the calculations, we have used the FLUENT 6.1 source code, which employs the volume of fluid (VOF) algorithm and a piecewise linear regularization (biviscosity model). We report quasistationary states that are established after relatively long transient periods. Results include the established velocity, interfacial shape, and extent of the yielded domain that depend strongly on the yield stress. Cases of coalescence and breakup of the drops are demonstrated. The yielded region enveloping the drops may include subdomains of floating unyielded material. Their location and shapes are very different from those encountered in the studies of solid inclusions in yield stress fluids. In section 2, the mathematical formulation of the problem is presented. The governing equation and boundary conditions describing the motion of a Newtonian drop in a tube, filled with a viscoplastic medium, are given. Difficulties encountered in numerical simulations of yield stress flow and the methods of regularization described in the literature are briefly discussed. Afterward, the method employed by FLUENT is described in more detail. Benchmark computations of the motion of a solid sphere are performed, and comparison with the results available in the literature is presented. Section 3 is devoted to the results of the performed simulations. In subsection 3.1, the motion of a single drop under the action of the gravity force inside a tube filled by a Bingham fluid is studied. Our computations revealed that, after a transient period, a quasistationary state in which the drop moves with constant velocity inside a liquid envelope of stationary shape is established. This established velocity, shape of the drop and yielded regions, and the time evolution of the velocity were chosen as characteristics of the flow. In subsection 3.1.1, we report the dependence of the sedimentation velocity on the auxiliary parameters of the numerical algorithm, such as the computational grid and the magnitude of the “viscosity of solid”. Subsection 3.1.2 is devoted to the description of the dependence of the established characteristics on the physical parameters of the problem. In subsection 3.2, the interaction of two drops was studied. Here, the process proved to be substantially nonstationary and our focus was on the dynamics of the deformation and coalescence of the drops. The numerical parameters (e.g., the mesh density and “solid state viscosity”) determined in the one-drop case were used for simulating the two-drop case as well. In section 4, we discuss the obtained physical results and the limits of the available software in the simulation of yield stress flows.
Figure 1. Schematic representation of the system of a drop falling in a tube filled with a viscoplastic medium.
where τ and D are the deviatoric stress and rate of strain tensors, while τY and ηp denote the yield stress and the “plastic viscosity”, respectively. |D| and |τ| are the second invariants of D and τ. The yield stress material with constitutive eqs 1 and 2 is addressed as a Bingham fluid. A popular generalization of such media is a Hershel-Bulkley fluid, with
(
τ ) K|D|(n-1) +
D ) 0, |τ| < τY
2.1. Modeling of Yield Stress Fluids. Yield stress, or viscoplastic, media deform as non-Newtonian viscous fluids when the local stress exceeds some critical value, while at lower levels of stress the material behaves as a solid. A properly invariant model for the deviatoric stress tensor formulated by Oldroyd10,11 and Prager12 employs a von Mises yield criterion for flow
(
(3) (4)
where n is a power law index and K is a consistency factor. When the power law index is unity, the Hershel-Bulkley model reduces to the Bingham model, while the consistency index is equivalent to the plastic viscosity. 2.2. Basic Equations and Boundary Conditions. Consider the gravity-induced motion of single or multiple viscous drops in a closed cylindrical tube of radius B and length L. The drops are placed on the axis of the tube and then fall under the action of gravity. When a particle moves within an otherwise quiescent viscoplastic medium, the induced stresses result in the appearance of a yielded region around the inclusion, while the rest of the ambient medium remains solid. The geometry of the flow induced by a drop falling inside a closed tube filled with a viscoplastic medium is schematically depicted in Figure 1. The flow on both sides of the drop’s surface at x ∈ Ωd and x ∈ Ωp, is governed by the conservation equations of mass and momentum
F 2. Mathematical Formulation of the Problem
)
τY D, |τ| g τY |D|
dV ) -∇p + F + ∇‚τ dt ∇‚V ) 0
(5) (6)
where V is the velocity vector, p is the scalar pressure, F is the density of the medium, and F is the density of the external body forces. In this paper, we consider gravity-induced motion, F ) Fg, where g is the gravitational acceleration. For a Newtonian drop, τ is proportional to the rate of deformation tensor,τ ) ηdD, where ηd is a constant viscosity. No-slip conditions are imposed on the tube boundary.
)
(1)
V)0
D ) 0, |τ| < τY
(2)
At the interface of the drop, conventional conditions of stress balance (dynamic boundary conditions)
τ ) ηp +
τY D, |τ| > τY |D|
(7)
6988
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006
[-pI + τ]‚n ) 2σHn
(8)
and continuity of velocity (kinematic boundary conditions)
[V] ) 0, V‚n ) Vn
(9)
are held. Here, [ ] denotes a difference across the surface, I is the unit tensor, σ is the surface tension coefficient, n is a unit normal vector, Vn is the velocity of the interface in the normal direction, and H is the mean curvature. The interfacial conditions (8)-(9) are valid in the case when the drop contacts the yielded region of the Bingham medium. For contact with the unyielded region, they should be substituted by the no-slip condition. The formulation of the problem is completed by specifying initial conditions for the shape and position of the drop and the velocity field. 2.3. Dimensionless Variables and Governing Parameters. At the outset, we introduce dimensionless variables based on the size of the smaller drop: a, for the length; V* ) a2∆F|g|/ ηp, for the velocity where ∆F ) Fd - Fp, with Fd and Fp being the density of the drop and of the ambient plastic medium, respectively; τ* ) ηpV*/a, for stress and pressure; and t* ) V*/a, for the time. The system is governed by the following set of dimensionless parameters: Reynolds number, Re ) FpaV*/ ηp; Bingham number, Bn ) τYa/(V*ηp); capillary number, Ca ) ηpV*/σ; density ratio, F ) Fd/Fp; viscosity ratio, λ ) ηd/ηp; and the geometric parameters, dimensionless radius, b ) B/a, and dimensionless length, l ) L/a, of the tube. In the case of multiple drops, their scaled dimensions, Ri ) ai/a, should be added to the set of governing parameters. In dimensionless variables, the governing equations take the form
Re
Fi dV ) -∇p + e + ∇‚τ, ∇‚V ) 0 dt ∆F x ∈ Ωi, i ) d or p
successfully used recently, e.g., by Vola et al.47 and by Frigaard and co-workers20,45 to model a number of single- and multiphase duct flows. Implementation of this approach to more general multiphase flows, especially to unsteady flows with surface tension that is the subject of our interest, is expected to encounter difficulties both in mathematics and numerics that are not yet resolved. The most popular alternative to the augmented Lagrangian approach is the use of regularized constitutive equations, in which the discontinuous Bingham model is replaced by a continuous one and where the solid is modeled by a liquid with a very high viscosity. One of the main advantages of this approach is the ease of its application in standard and commercial solvers for partial differential equations. The regularized models differ in the adopted rheology of the solid-state liquid. The most popular approaches, as in the work of Beris et al.32 and Liu et al.,31,40 used a regularization parameter, proposed by Bercovier and Engelman48
(
τ) 1+
)
Bn D |D| + PR
(12)
A different regularization model is the exponential expression
(
τ) 1+
Bn [1 - exp(-m|D|)] D |D|
)
(13)
which was employed by, e.g., Papanastasiou,49 Mitsoulis and co-workers,29,30 Alexandrou et al.50 and Li and Renardy.46 These models, characterized by regularization parameters, PR and m, approximate the ideal Bingham fluid constitutive eqs 1 and 2 at PR f 0 and m f ∞. Yet another approximation is the biviscosity model with a piecewise linear constitutive equation14 that is now commonly used, e.g., in the software FLUENT, where
(10)
τ ) ηeffD
(14)
where e ) g/|g|, while
(
ηeff ) 1 +
Bn D, |τ| > Bn, x ∈ Ωp τ) 1+ |D|
(
)
[
])
ηeff ) η, |D| e
D ) 0, |τ| < Bn, x ∈ Ωp τ ) λD, x ∈ Ωd
Bn 1 Bn 1 - , |D| > |D| η η
(11)
The dynamic boundary conditions at the interface of the drop become
Ca[-pI + τ]‚n ) 2Hn while the kinematic conditions at the interface, the no-slip conditions at the tubes surface, and the initial conditions retain their original form. 2.4. Numerical Simulations of Bingham Fluid Flow. Discontinuity in the constitutive equation, separating the unyielded and yielded states, creates significant difficulties in the numerical solution of viscoplastic problems. The two main general approaches to the numerical modeling of such media are the augmented Lagrangian approach and viscosity regularization. Several hybrid methods have been developed, which proved successful in capturing both yielded and unyielded regions for specific problems. The augmented Lagrangian approach is based on variational formulation of the problem,18,21 and it allows one to cope with the exact discontinuous constitutive law. This approach was
Bn η
(15) (16)
For low strain rates (|D| e Bn/η), the “rigid” material acts like a very viscous Newtonian fluid with viscosity η that plays the role of a regularization parameter. Simulations performed for various values of η revealed a convergence of the results for the velocity field and the form of the yielded domain at η f ∞. The results discussed below were obtained at a fixed value of η, the choice of which is discussed in section 3. In what follows, the flow region where (16) is satisfied is addressed as an unyielded region, though the velocity field does not actually vanish there. The common practice is that the use of regularization facilitates the calculation and convergence of the velocity field but provides a rather inaccurate prediction of the stress field. Indeed, a recent study by Frigaard and Nouar51 compared the rate of convergence of various two-dimensional stress fields using several regularization rheological models and found the one suggested by Papanastasiou49 to be superior and the biviscosity models to lag behind. Thus, the results below concerning the exact shapes of the yielded domain and the yield limit should be viewed with some caution and are awaiting further corroboration by other methods of calculation.
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006 6989 Table 1. Comparison of Drug Coefficient Results with Those Computed Making Use of the Empirical Formula (16) from Ref 30 Bn
CS computed with Fluent
CS computed with empirical formula from ref 30
0 0.108 0.747 2.299 8.047 59.59 544.6
1.946 2.19 3.52 6.35 15.11 82.61 698.8
1.98 2.23 3.45 6.11 15.06 84.48 633.8
It should be noted also that the definition of Reynolds number using ηp is rather arbitrary and does not reflect fully the nature of the flow, since the actual viscosity in the yielded region varies considerably between ηp and the regularization rigid material viscosity η. Hence, while the calculated values of Re appear relatively high (the results presented in section 3 were computed at Re ) 3), the motion is dominated by viscous stresses. Our simulations were performed using the source code FLUENT 6.1 which employs a volume of fluid (VOF) numerical scheme. Use was made of a uniform mesh with square cells. Most of the presented results were obtained employing a mesh with 64 000 cells. Note that the VOF method has an inherent inaccuracy since it reconstructs the interface after computing the velocity and pressure fields considering the drop and its surroundings as a single phase. This makes the interface rather diffusive in each computational cell. As a benchmark, we chose the problem of a solid sphere falling in a tube filled with a Bingham fluid30 since no results on the motion of deformable drops in a partially yielded Bingham fluid are available so far, and the crucial source of the computational difficulties in the problem under consideration is a highly nonlinear and discontinuous constitutive equation for the Bingham fluid. Note that the results of ref 30 were recently replicated in ref 31, making use of another numerical method. Following refs 30 and 31, we calculated stationary flow around a spherical particle, which moves along the center of a cylindrical tube with a given constant velocity, V. A reference frame linked to the particle was chosen. The velocity was scaled by V. Thus, the velocity was set zero at the boundary of the sphere, while at the top, bottom, and side walls the velocity is directed along the axis, having a unit absolute value. The total force, F, exerted on the particle by the flow was computed with the Stokes drag coefficient defined as
CS )
F 6πµVR
The drag coefficients of ref 30 are reproduced fairly well in our computations as is illustrated in Table 1, where the drag coefficient is computed for the radii ratio R ) 4 (this R is used further in our original computations with the falling drop). A similar reproducibility was obtained for the yield surfaces. 3. Results of Simulations Our calculations were performed for the gravity-induced motion in a closed cylindrical tube of radius b ) 4 and of length l ) 40. One or two initially spherical drops of a density F ) 1.54 and viscosity λ ) 0.07 are placed on the axis of symmetry of the tube and then fall under the action of gravity. The tube is filled with a Bingham fluid. All the results presented below are computed for Re ) 3, while the values of the Bingham and capillary numbers vary. Our computations were performed assuming an axial symmetry. The stability of such configuration was not studied. It is
Figure 2. Temporal evolution of the velocity of the center of mass of a single drop falling in Bingham fluid. Ca ) 17.33. The Bingham number in the range 0.058-0.13.
Figure 3. Temporal evolution of the velocity of the center of mass of a single drop falling in Bingham fluid. Bn ) 0.058. The capillary number in the range 1.7-8.7.
well-known that, for a somewhat related problem of a pairwise interaction of drops in an unbounded Newtonian fluid, an axisymmetric configuration with the larger drop following the small one is unstable and even minor disturbances result in overtaking of the small particle by the larger one. Thus, the results of axisymmetric computations should be regarded with caution especially when the interaction of the two drops is concerned. On the other hand, the presence of the tube wall confining the flow region and the fact that the drops in our computations are of equal size are expected to have a stabilizing effect on the axisymmetric configuration. 3.1. Single-Drop Motion. Results of the simulation of singledrop gravity-induced motion are illustrated in Figures 2-10. Initially, the drop is separated by a distance of two diameters from the upper solid edge of the pipe. The temporal evolution of the velocity of the falling drop’s center of mass is presented in Figures 2 and 3. In Figure 2, Ca ) 17.33 and various curves correspond to various values of the Bingham number, while, in Figure 3, Bn ) 0.058 and various curves correspond to various values of the capillary number. One can see that in all cases after a transient period of increase the velocity reaches some maximum, then decreases, and tends to an almost stationary value. The drop and the yielded region assume an established shape after the same period of time. However, a careful examination in a refined scale reveals that a steady state is not achieved in a rigorous sense but that the velocity is constantly decreasing (see Figure 3). We attribute this affect to a finite length of a tube in our computations. It is anticipated that these unsteady effects will decay with the increase of the pipe’s length. 3.1.1. Study of the Convergence. We chose the time evolution of the drop’s center of mass, its established velocity,
6990
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006
Figure 4. Established velocity of the drop versus solid state viscosity η, for Ca ) 17.33. The upper, middle, and lower curves are calculated for Bn ) 0.077, 0.096, and 0.12, respectively.
Figure 5. Established velocity of a single drop falling in Bingham fluid versus the number of computational elements. Ca ) 17.33. The upper, middle, and lower curves are calculated for Bn ) 0.077, 0.096, and 0.12, respectively.
and the established shape of the drop to characterize the motion and to check the convergence of our computations with respect to the regularization parameter of the model and the refinement of the computational grid. This check was performed for several values of the input parameters. The results of the checks are illustrated below for the established velocity, and its dependence on the value of regularization parameter η and the computational mesh density is shown in Figures 4 and 5. In Figure 4, the established velocity is plotted versus the solid state viscosity η for Ca ) 17.33. The upper, middle, and lower curves are calculated for Bn ) 0.077, 0.096, and 0.12, respectively. Convergence of the results is evident at low Bn cases (two upper curves), where the velocity does not change for η > 3400, while, for higher Bn, the convergence is considerably slower and the limit, probably, is not fully achieved for η < 7000. The effect of η on the shapes of the drops and the yielded regions is discussed in section 3.1.2. All the computations presented henceforth are performed for η ) 3400. The convergence of the method with mesh refinement is illustrated in Figure 5, where the established velocity is plotted versus the number of cells in the computational grid, N, for Ca ) 17.33. The upper, middle, and lower curves are calculated for Bn ) 0.077, 0.096, and 0.12, respectively. A tendency to convergence at high N is evident at every Bn. The results presented below were obtained with N ) 64 000. These results are expected to be close to the limiting values and, on the other hand, are obtained in a reasonable computational time. 3.1.2. Effect of Physical Parameters. In this subsection, we describe the influence of the physical parameters on the process. The computational grid and regularization parameter η are the same in all the runs. The duration of the transient period as well as the established velocity depends strongly on the Bingham number, while the dependence on the capillary number is
Figure 6. Established velocity of a single drop falling in Bingham fluid versus capillary number at various Bingham numbers in the range 0.0580.12.
Figure 7. Established velocity of a single drop falling in Bingham fluid as a function of the Bingham number at Ca ) 17.33.
considerably weaker. The latter dependence is illustrated in Figure 6. It is evident that the velocity grows with the capillary number starting with some finite positive value at Ca ) 0, corresponding to the case of nondeformable drop, and tends to another finite value at Ca f ∞, corresponding to the case of vanishing interfacial tension. The ratio between these two limiting numbers is ∼1.2, reflecting the observation that interfacial deformations are moderate. The dependence of the established velocity on the Bingham number for Ca ) 17.33 is illustrated in Figure 7. It is evident that the velocity decays fast with the growth of Bn and vanishes when it reaches some finite critical value. This effect can be explained as follows: The established velocity is determined by the balance of gravity and viscous resistance forces. The latter depends strongly on the flow domain that, in turn, is a function of the yield stress. The larger Bn is the smaller the yielded domain is and the larger the resistance is; hence, the smaller the velocity is. It can be anticipated that some critical value of Bn exists, above which the induced stresses are not strong enough to yield the Bingham medium adjacent to the drop surface and, hence, in such cases, no flow occurs. Note again that in the performed calculations zero flow actually never occurs. For high Bn, the anticipated characterictic velocities are determined by the Newtonian flow in a “solid” region, being on the order of 1/η. Thus, the computed critical value of Bn is probably overestimated and provides only a crude approximation to the yield limit. Figure 8 illustrates the dependence of the form and dimension of the yielded region on the magnitude of yield stress. The local viscosity of the media is marked by colors. The very low viscosity blue domain corresponds to the Newtonian drop; the unyielded solid is marked by red, while the colored region between these two domains corresponds to yielded Bingham material with variable viscosity.
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006 6991
Figure 8. Samples of the established shapes of a single drop and the yielded domain for various Bingham and capillary numbers: (left column) Ca ) 52, (right column) Ca ) 2.9; (upper row) Bn ) 0.058, (middle row) Bn ) 0.12, (lower row) Bn ) 0.15. The value of the right-hand side (RHS) scale is log ηeff.
Figure 9. Streamlines of established flow induced by a single drop falling in a Bingham fluid. Ca ) 52, and Bn ) 0.058. The color scale is the same as that in Figure 8.
For relatively small yield stress (see the first row of Figure 8), the drop moves in a wide liquid domain and its shape is similar to that of a drop in a Newtonian fluid. In the first row of Figure 8, one can distinguish a peach-shaped domain with relatively low effective viscosity, a vicinity of the boundary of this domain with the viscosity changing fast in the normal direction, small fuzzy regions in the yielded domain with large variable viscosity up and downstream of the peach domain, and an unyielded region with a constant viscosity equal to that of the solid state. A closer look at Figure 8a is given in Figure 9, where the streamlines in the laboratory reference frame are presented as well, indicating downward motion of the drop. The streamline pattern differs considerably from that induced by the drop falling in a tube filled with a Newtonian fluid. In the latter case, the streamlines are open, nearly straight lines, while for a Bingham fluid the streamlines are closed, being probably similar to the streamlines induced by the drop falling inside the cavity shaped like the yielded region.
Figure 10. Components of the scaled strain rate tensor. Dzz, Drr, and Drz are presented in the upper, middle, and lower row, respectively. Ca ) 52. (a) Bn ) 0.058. (b) Bn ) 0.12.
The toroidal vortex that is seen in Figure 9 assumes a stagnation ring in its center. In the vicinity of this stagnation ring, the rate of strain is lower and, hence, the effective viscosity is higher than in the neighboring regions. These domains of high viscosity are almost unseen in Figure 8 but can be distinguished in Figure 9. In the second row of Figure 8, where the Bingham number is higher, the peach-shaped domain shrinks, while the fuzzy regions grow. The toroidal domain of high viscosity in the center of the vortex is more pronounced. The yielded domain is considerably smaller, and the effective viscosity in this region is higher; hence, the hydrodynamic resistance to the translation of the drop is much larger than in the previous case and the established velocity is lower. For high Bn (last row in Figure 8), a solid zone of unyielded material comes close to the drop sides, squeezes the drop, and results in a considerable deceleration of the motion. It is evident that the fuzzy regions, in which the variable viscosity is large and close to the solid state viscosity, are considerably extended, while the low viscosity domain shrinks to the vicinity of the drop’s surface. Thus, the drop moves in a small envelope of yielded fluid with high viscosity and, as a result, the sedimentation velocity becomes vanishingly small. For larger yield stress, the Bingham fluid is hardly yielded, corresponding to the cases of vanishing velocity shown in Figure 7. The components of the scaled strain rate tensor, D, are shown in Figure 10. The left and right columns depict the components associated with the drops in Figure 8a and c, respectively. The top and middle rows show the two normal diagonal components, Dzz and Drr, while the bottom row shows the axisymmetric shear component, Drz. The relative nondimensional intensity is given by the color scale. Note the compression and extension regions near the front and rear poles and the strong shear component around the drop equator. Note also the intensification of the equatorial shear components with the increase of the Bingham number and the shrinking diameter of the yielded region. In the left-hand side, one can clearly detect the contours showing
6992
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006
Figure 11. Temporal evolution a of pair of highly deformable drops falling in Bingham fluid. Bn ) 0.058, and Ca ) 52.
the steep changes in velocity gradients in the vicinity of the transition from the yielded to unyielded regions. It is also clear that, in view of the regularization of the rheological model, some velocity variations can be detected beyond that contour within the unyielded domain. 3.2. Interaction of Two Drops. Two equal spherical drops, R ) 1, are initially placed at the axis of the cylindrical tube (the upper drop’s center is located four radii apart from the upper edge of the pipe and the center of the lower one is four radii apart from the upper drop center). The drops move under the action of gravity and buoyancy forces. A series of simulations for various values of Bingham and capillary numbers have been performed. Samples of the pairwise interaction of the drops are demonstrated in Figures 11-14. In Figure 11, Bn ) 0.058 and Ca ) 52; in Figure 12, Bn ) 0.12 and Ca ) 52; in Figure 13, Bn ) 0.058 and Ca ) 1; in Figure 14, Bn ) 0.12 and Ca ) 1. The various colors correspond to various values of the local effective viscosity, as in Figure 8. Drops are dark blue, while the unyielded domain is red. Strong variations of local viscosity are evident inside the yielded region. The form and dimension of the yielded region is determined by the gravity force acting on the drops and the magnitude of the yield stress, while the form of the drops is a result of the interplay of viscous and capillary forces. Initially, the configuration of the two drops, the form of the yielded domain, and the shape of the drops are almost symmetric with respect to a plane normal to the centerline in its center. A
Figure 12. Temporal evolution of a pair of highly deformable drops falling in Bingham fluid. Bn ) 0.12, and Ca ) 52.
toroidal domain of high viscosity and a thin island of unyielded material is seen in the vicinity of the center plane. With the passage of time, the drops move and deform with the leading and the trailing drops attaining oblate and prolate shapes, respectively. Thus, the hydrodynamic resistance to the motion of the leading drop increases and its motion slows down, while the resistance to the motion of the trailing one decreases causing advance toward the leading drop. As a result, the drops approach each other and eventually merge into one large drop. As the drops approach each other, the domain of high viscosity shrinks and the solid island “melts”. The details of the approach process and the further behavior of the large drop depend strongly on the Bn and Ca numbers. For highly deformable drops with Ca ) 52, corresponding to almost zero surface tension (Figures 11 and 12), the leading drop attains a crescent shape with sharp edges, while the trailing one elongates considerably forming a long tail with a sharp end. With the further passage of time, the trailing drop penetrates the leading one together with the portion of Bingham fluid in front of it. When the two drops merge, this portion of Bingham fluid is found encapsulated inside the newly formed large drop, while the tail elongates and breaks up forming a set of daughter drops. It should be noted that since the tracking of the interface, using VOF method, tends to describe the surfaces in a somewhat diffuse manner, dramatic changes such as coalescence, breakup, and tearing should be viewed with caution. For lower Ca (Figures 13 and 14), the drops are less deformable and attain a rounded shape up to the moment of coalescence, since the considerable surface tension does not
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006 6993
Figure 13. Temporal evolution of a pair of drops falling in Bingham fluid. Bn ) 0.058, and Ca ) 1.
Figure 14. Temporal evolution of a pair of drops falling in Bingham fluid. Bn ) 0.12, and Ca ) 1.
allow the formation of sharp angles and tails. After the period of approach, the drops evidently touch at one centerline point and coalesce forming one large drop that, in this case, does not contain any inclusion. No breakup or formation of daughter drops take place, and after a relatively large transient period, the drop attains almost a spherical stationary shape. Variation of the Bingham number between 0.05 and 0.2 does not change the scenario of drop interaction, but its effect on the velocities is quite strong. For high Bn (Figures 12 and 14), the region of low viscosity in the yielded domain shrinks to the relatively small vicinity of the drops and the toroidal solid island is considerably extended. As a result, the motion is much slower. Note that there is a difference in the sedimentation speed between the pairs that are shown in Figures 11 and 13 (for the same Bingham number). This difference may be attributed to two effects. There is a small effect of the different drop deformation patterns for the different Ca values on the hydrodynamic resistance to the motion. This effect is operative in the single-drop case as well (see Figure 6). However, the major effect is the different settling velocities of a coalesced bigger drop and of two drops settling in close proximity. In Figure 15, the time in which the two drops come into contact is plotted versus ln(Ca) for various values of Bn. An extremely strong dependence on both governing parameters is evident at small Ca, while, at high Ca, the approach seems to exhibit a small nonzero limiting value, which depends slightly on Bn.
Figure 15. Time in which the two drops come into contact versus ln(Ca) for various values of the Bingham number Bn.
4. Conclusions Numerical simulation of the gravity-induced motion of single and multiple deformable Newtonian drops in tubes filled with Bingham plastic were performed using the FLUENT 6.1 source code, which employs the volume of fluid (VOF) algorithm and piecewise linear regularization of the constitutive equation. One or two drops were initially placed at the axis of a cylindrical tube and moved under the action of gravity and buoyancy forces. The focus was on the influence of the drop deformations on the shape and extent of the yield regions and on the velocities of the drops.
6994
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006
For a single drop, a quasistationary state is established after a relatively long transient period. The duration of the transient period and the established velocity depend strongly on the Bingham number. At higher Bn, a smaller domain is yielded, that leads to a larger viscous resistance and, hence, a smaller velocity. There is some critical value of Bn, above which the induced stresses are not strong enough to yield the Bingham medium and, hence, no flow occurs. The dependence on the capillary number is less pronounced. The interaction of two drops results in a substantial extension of the yielded domain and in the acceleration of the motion. Two equal-size drops tend to approach each other and coalesce. The interaction pattern depends strongly on the deformability of the drops, which is characterized by the capillary number. Higher deformations result in higher stresses and, thus, in an extended yielded region and higher migration velocities. After coalescence, a single large fluid particle is formed that at high Ca breaks up forming one large leading drop and one or several small trailing ones. For low Ca, the newly formed drop does not break up but attains, after a transient period, a stationary, nearly spherical shape. Acknowledgment This work was partially supported by the Asher Space Research Center at the Technion. O.M.L. and A.P. acknowledge the support of the Center for Absorption in Science, Ministry of Immigrant Absorption State of Israel. Nomenclature a ) radius of the drop B ) radius of the tube b ) scaled radius of the tube Bn ) Bingham number Ca ) capillary number D ) rate of strain tensor |D| ) the second invariant of D d ) initial separation distance F ) density of the external body force H ) the mean curvature I ) unit tensor K ) consistency factor L ) the length of the tube l ) scaled length of the tube m ) regularization parameter n ) unit normal vector n ) power law index p ) pressure PR ) regularization parameter Re ) Reynolds number V ) velocity vector Vn ) velocity of the interface in the normal direction Rq ) volume fraction η ) constant viscosity of “rigid state” ηd ) constant viscosity of the drop ηeff ) variable effective viscosity of Bingham fluid ηp ) the plastic viscosity σ ) interfacial tension λ ) viscosity ration F ) density ratio Fp ) density of the medium Fd ) density of the drop τ ) deviatoric stress tensor
|τ| ) the second invariant of τ τY ) the yield stress Literature Cited (1) Bird, R. B.; Dai, G. C.; Yarusso, B. J. The rheology and flow of viscoplastic materials. ReV. Chem. Eng. 1983, 1, 1-70. (2) Barnes, H. A. The yield stress - a review or ′πRντR Fι′ everything flows? J. Non-Newtonian Fluid Mech. 1999, 81, 133-178. (3) Coussot, P. Mudflow rheology and dynamics; Balkema: Rotterdam, 1997. (4) Nielloud, F.; Mart-Mestres, G. Pharmaceutical emulsions and suspensions; Marcel Dekker: New York, 2000. (5) Johnson, S. B.; Franks, G. V.; Scales, P. J.; Boger, D. V.; Healy, T. W. Surface chemistry-rheology relationships in concentrated mineral suspensions. Int. J. Mineral Process. 2000, 58, 267-304. (6) Zhou, Z.; Solomon, M. J.; Scales, P. J.; Boger, D. V. The yield stress of concentrated flocculated suspensions of size distributed particles. J. Rheol. 1999, 43, 651-671. (7) Kang, J.; Butler, P. B.; Baer, M. R. A thermodynamical analysis of hot spot formation in condensed-phase, energetic materials. Combust. Flame 1992, 89, 117-139. (8) Gallardo, V.; Ruiz, M. A.; Delgado, A. V. Pharmaceutical suspensions and their applications. In Pharmaceutical emulsions and suspensions; Nielloud, F., Mart-Mestres, G., Eds.; Marcel Dekker: New York, 2000. (9) Bingham, E. C. Fluidity and Plasticity; McGraw-Hill: New York, 1922. (10) Oldroyd, J. D. A rational formulation of the equations of plastic flow for a Bingham solid. Proc. Cambridge Philos. Soc. 1947, 43, 100105. (11) Oldroyd, J. D. Two-dimensional plastic flow of a Bingham solid. Proc. Cambridge Philos. Soc. 1947, 43, 383-395. (12) Prager, W. On slow viscoplastic flow. In Studies in Mathematics and Mechanics: R. Von Mises presentation Volume; Academic: New York, 1954; pp 208-216. (13) Burgos, G. R.; Alexandrou, A. N.; Entov, V. On determination of yield surfaces in Hershel-Bulkley fluids. J. Non-Newtonian Fluid Mech. 1999, 43, 463-483. (14) Lipscomb, G. G.; Denn, M. M. Flow of Bingham fluids in complex geometries. J. Non-Newtonian Fluid Mech. 1984, 14, 337-346. (15) Mosolov, P. P.; Miasnikov, V. P. Variational methods in the theory of the fluidity of a viscous-plastic medium. J. Mech. Appl. Math. 1965, 29, 468-492. (16) Yoshioka, N.; Adachi, K. On variational principles for a nonNewtonian fluid. J. Chem. Eng. Jpn. 1971, 4, 217-220. (17) Duvaut, G.; Lions, J. L. Inequalitie in Mechanics and Physics; Springer-Verlag: New York, 1976. (18) Glowinski, R. Numerical Methods for Nonlinear Variational Problems; Springer series in computational physics; Springer-Verlag: New York, 1984. (19) Huilgol, R. R. Variational principles and variational inequalities for a yield stress fluid in the presence of slip. J. Non-Newtonian Fluid Mech. 1998, 75, 231-251. (20) Frigaard, I. A.; Leimgruber, S.; Scherzer, O. Variatinal methods and maximum residual wall layers. J. Fluid Mech. 2003, 483, 37-65. (21) Glowinski, R.; Lions, J. L.; Tremolieres, R. Numerical Analysis of Variational Inequalities; North-Holland Publishing Company: Amsterdam, 1981; Vol. 8. (22) Nguyen, Q. D.; Boger, D. V. Yield stress measurements for concentrated suspensions. J. Rheol. 1983, 27, 321-336. (23) Nguyen, Q. D.; Boger, D. V. Direct yield stress measurement with the vane method. J. Rheol. 1985, 29, 335-347. (24) Spaans, R. D.; Williams, M. C. Letter to the Editor: At last: a true liquid-phase yield stress. J. Rheol. 1995, 39, 241-246. (25) Santos, O. L. A.; Azar, J. J. A study of gas migration in stagnant non-Newtonian fluids; Technical Report SPE 39019, Society of Petroleum Engineers: Richardson, TX, 1997. (26) Andres, V. T. Equilibrium and motion of a sphere in a visco-plastic fluid. Dokl. Acad. Nauk SSSR 1960, 133, 777-780. (27) Valentic, L.; Whitmore, R. L. The terminal velocity of spheres in Bingham plastics. Brit. J. Appl. Phys. 1965, 16, 1197-1203. (28) Yoshioka, N.; Adachi, K.; Ishimura H. On creeping motion of viscoplastic fluid past a sphere. Kagaku Kogaku 1971, 10, 1144-1152. (29) Blackery, J.; Mitsoulis, E. Creeping motion of a sphere in tubes filled with a Bingham plastic. J. Non-Newtonian Fluid Mech. 1997, 70, 59-77.
Ind. Eng. Chem. Res., Vol. 45, No. 21, 2006 6995 (30) Beauline, M.; Mitsoulis, E. Creeping motion of a sphere in tubes filled with Herschel-Bulkley fluids. J. Non-Newtonian Fluid Mech. 1997, 72, 55-71. (31) Liu, B. T.; Muller, S. J.; Denn, M. M. Convergence of a regularization method for creeping flow of a Bingham material about a rigid sphere. J. Non-Newtonian Fluid Mech. 2002, 102, 179-191. (32) Beris, A. N.; Tsamopoulos, J. A.; Armstrong, R. C.; Brown R. A. Creeping motion of a sphere through a Bingham plastic. J. Fluid Mech. 1985, 158, 219-244. (33) Roquet, N.; Saramito, P. An adaptive finite element method for Bingham fluid flow around a cylinder. Comput. Methods Appl. Mech. Eng. 2003, 192, 3317-3341. (34) Smyrnaios, D. N.; Tsamopoulos, J. A. Squeeze flow of Bingham plastics. J. Non-Newtonian Fluid Mech. 2001, 100, 165-190. (35) Ansley, R. W.; Smith, T. N. Motion of a spherical particle in a Bingham plastic. AIChE J. 1967, 13, 1193-1196. (36) de Plessis, M. P.; Ansley, R. W. Settling the parameters in solid pipelines. J. Pipeline DiV. ASCE 1967, 93, 1-17. (37) Ito, S.; Kajachi, T. Drag force on a sphere moving in a plastic solid. J. Chem. Eng. Jpn. 1969, 2, 19-24. (38) Atapattu, D. D.; Chabra, B. P.; Uhlherr, P. H. T. Creeping sphere motion in Herschel-Bulkley fluids. Flow field and drag. J. Non-Newtonian Fluid Mech. 1995, 59, 245-265. (39) Jossic, L.; Magnin, A. Drag and stability of objects in a yield stress fluid. AIChE J. 2001, 47, 2666-2672. (40) Liu, B. T.; Muller, S. J.; Denn, M. M. Interactions of two rigid spheres translating collinearly in creeping flow in a Bingham material, J. Non-Newtonian Fluid Mech. 2003, 113, 49-67. (41) Jossic, L.; Magnin, A. Structuring under flow of suspensions in a gel. AIChE J. 2004, 50, 2691-2696.
(42) Jossic, L.; Magnin, A. Structuring of gelled suspensions flowing through a sudden three-dimensional expansion. J. Non-Newtonian Fluid Mech. 2005, 127, 201-212. (43) Bhavaraju, S. M.; Mashelkar, R. A.; Blanch, H. W. Bubble motion and mass transfer in Non-Newtonian fluids. AIChE J. 1978, 24, 10631070. (44) Stein, S.; Buggisch, H. Rise of pulsating bubbles in fluids with a yield stress. Z. Angew. Math. Mech. 2000, 80, 827-834. (45) Dubash, N.; Frigaard, I. Conditions for static bubbles in viscoplastic fluid. Phys. Fluids 2004, 16, 4319-4330. (46) Li, J.; Renardy, Y. Y. Shear-induced rupturing of a viscous drop in Bingham liquid. J. Non-Newtonian Fluid Mech. 2000, 95, 235-251. (47) Vola, D.; Boscardin, L.; Latche, J. Laminar unsteady flows of Bingham fluids: A numerical strategy and some benchmark results. J. Comput. Phys. 2003, 187, 441-456. (48) Bercovier, M.; Engelman, M. A finite element method for incompressible non-Newtonian flows. J. Comput. Phys. 1980, 36, 315-326. (49) Papanastasiou, T. C. Flow of materials with yield. J. Rheol. 1987, 31, 385-404. (50) Alexandrou, A. N.; McGilvreay, T. M.; Burgos, G. R. Steady Herschel-Bulkley fluid in tree dimensional expansions. J. Non-Newtonian Fluid Mech. 2001, 100, 77-96. (51) Frigaard, I.; Nouar, C. On the usage of viscosity regularization methods for visco-plastic fluid flow. J. Non-Newtonian Fluid Mech. 2005, 127, 1-26.
ReceiVed for reView November 3, 2005 ReVised manuscript receiVed January 26, 2006 Accepted March 20, 2006 IE051222E
