Transport and dispersion ofparticles in visco-plastic fluidsbySarah HormoziB.Sc. Mechanical Engineering, Shiraz University, 2003M.Sc. Mechanical Engineering, Sharif University of Technology, 2006Ph.D. Mechanical Engineering, The University of British Columbia, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Mathematics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)August 2013c? Sarah Hormozi 2013AbstractThis thesis focuses on development of a model to predict ?spreading? of thesolids (i.e. proppant) fraction during the fracturing operation. We developa 1D model that allows us to estimate dispersion of solid particles along avertical pipe in a fully turbulent flow of a shear thinning yield stress fluid (i.e.,visco-plastic fluid), as well as slip relative to the mean flow. In dimensionlessform, this results in a quasilinear advection-diffusion equation. Advection bythe mean flow, particle settling relative to the mean, in the direction of gravity,turbulent particle dispersivity and Taylor dispersion are the 4 main transportphenomena modelled in the 1D model. We provide a simple analysis of the1D model, suitable for spreadsheet-type field design purposes, in which weestimate ?mixing lengths? due to both settling and dispersion. Secondly, weprovide an accurate numerical algorithm for solution of the 1D model andshow how pulses of proppant (i.e. slugs) may or may not interact for typicalprocess parameters.iiPrefaceThe author of this thesis was the principal contributor to this research. Pro-fessor Ian Frigaard supervised the research.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Objectives of the thesis . . . . . . . . . . . . . . . . . . . . . 11.2 Literature overview . . . . . . . . . . . . . . . . . . . . . . . 31.3 Dimensional analysis . . . . . . . . . . . . . . . . . . . . . . . 51.4 Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . 92 Modeling approach . . . . . . . . . . . . . . . . . . . . . . . . . 102.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.1.1 Mass conservation and particle diffusivity . . . . . . . 102.1.2 Momentum balance . . . . . . . . . . . . . . . . . . . 122.1.3 Viscous stress closure . . . . . . . . . . . . . . . . . . 132.1.4 Turbulent velocity closure . . . . . . . . . . . . . . . . 16ivTable of Contents2.1.5 Relative velocity closure . . . . . . . . . . . . . . . . . 232.2 Particle dispersion . . . . . . . . . . . . . . . . . . . . . . . . 252.2.1 Dispersion model . . . . . . . . . . . . . . . . . . . . . 282.2.2 Particle diffusivity models . . . . . . . . . . . . . . . . 342.2.3 Slug spreading estimates . . . . . . . . . . . . . . . . . 403 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 433.1 Algorithm development . . . . . . . . . . . . . . . . . . . . . 433.1.1 Advection equation . . . . . . . . . . . . . . . . . . . 443.1.2 Diffusion equation . . . . . . . . . . . . . . . . . . . . 473.1.3 Advection-Diffusion equation . . . . . . . . . . . . . . 523.2 Physical results . . . . . . . . . . . . . . . . . . . . . . . . . . 524 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63vList of Tables1.1 Typical ranges of the dimensional process parameters in theindustrial application. . . . . . . . . . . . . . . . . . . . . . . 71.2 The range of non-dimensional parameters in the industrial ap-plication. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7viList of Figures1.1 Schematic of hydraulic fracturing . . . . . . . . . . . . . . . . 22.1 Amplification of the yield stress and consistency within the sus-pension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162.2 Dependency of bulk Re and B on mean flow and colormap offriction factor in (U?0, ??) plane . . . . . . . . . . . . . . . . . . 202.3 Dependency of bulk Re and B on power law index n and col-ormap of friction factor in (n, ??) plane . . . . . . . . . . . . . 212.4 Colourmaps of various flow parameters plotted in the (U?0, ??)plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 262.5 Colourmaps of various flow parameters plotted in the (dp, ??) plane 272.6 Variations in Pe for the different particle diffusivity models . . 382.7 Variations in DTPe2 for the different particle diffusivity models 392.8 Estimates of the lengths that a proppant slug will spread, basedon particle settling (slip), Lmix1, and particle diffusion, Lmix2 . 423.1 Solution of the linear (constant speed) advection problem in aperiodic domain with unit length by the NT scheme . . . . . . 483.2 Solution of the linear (constant speed) advection problem in aperiodic domain with unit length by the NT scheme . . . . . . 483.3 Numerical solution of nonlinear diffusion problem . . . . . . . 503.4 Numerical solution of nonlinear diffusion problem . . . . . . . 51viiList of Figures3.5 Exact (solid line) and numerical solution (dashed line) of thelinear (constant coefficient) advection-diffusion equation at dif-ferent times . . . . . . . . . . . . . . . . . . . . . . . . . . . . 533.6 The evolution of a square wave at different time under transportequation (2.70) . . . . . . . . . . . . . . . . . . . . . . . . . . 543.7 The interaction of two square waves transporting with equation(2.70) at different time . . . . . . . . . . . . . . . . . . . . . . 543.8 Results of solving (2.70) for a typical set of process parameters 563.9 A typical example on spreading of particles in the the pipe . . 573.10 A typical example on spreading of particles in the the pipe . . 583.11 A typical example on spreading of particles in the the pipe . . 59viiiAcknowledgmentsI would like to express my sincere gratitude and appreciation to those whoassisted me during the course of this program.Professor Ian Frigaard for serving as my advisor since 2008. Ian gave methe opportunity to work in several different areas of fluid mechanics rangingfrom hard core computations to lab-scale experiments to applied mathematicaltechniques. Ian showed me the beauty of Mathematics through his insightfulscientific guidance. This inspired me to persue my education in Mathematics tobridge science and engineering. Ian?s constant support and amiable personalityhave been of great value to me. I was truly blessed to be his student duringthe best school years of my life at UBC.Dr Dmitry Eskin for his useful advice on this research.Professor Neil Balmforth, Professor Michael Ward and Professor GeorgeBluman for their supervision, advice, and guidance on some other researchdirections. Additionally, I was fortunate to take few Applied Mathematicscourses from these gifted teachers.Professor Bud Homsy for his guidance and support throughout, for keepingan open door, for helpful advice on presentation and for shared physical insight.My friends in UBC and elsewhere, for their company and friendship. Inparticular, Dr Ali Vakil, Dr Mona Rahmani, Dr Anirban Guha, ProfessorBahman Gharesifard and Mahtab Fathinejad for their continuous support.My close family and my husband, Hamed, who stood by me and supportedall my ideas and dreams. I love you dearly.Finally, financial support of Schlumberger is gratefully acknowledged.ixDedicationTo my family and friends.xChapter 1IntroductionIt is common that many oil and gas wells do not produce at the expectedrate. Hydraulic fracturing is a stimulation technique used to increase theproductivity of wells in oil and gas reservoirs. In hydraulic fracturing, speciallyengineered suspensions are pumped at high pressure and rate into the reservoir,causing a propagating fracture to open. When the pressure is released thefracture is supported by the grains of solid (called proppant) that are leftbehind (see Fig.1.1). The increase in productivity results because the hydraulicconductivity of the fracture, filled with proppant, is significantly higher thanthat of the surrounding formation. More details on hydraulic fracturing canbe found in [46].A recent trend in the oil industry is to use the cyclic pumping of a prop-pant slurry interspersed with clear fracking fluid. This procedure is found toincrease the subsequent productivity (see e.g., [11, 39, 40]), over that due toconventional fracturing with a continuous stream of propant. It is thereforeof interest to understand how slugs of proppant pumped in a cyclic fashioncan disperse, both in the pipe on the way to the fracture and within the frac-ture itself. More clearly, under which physical situations do slugs of proppantinterfere with one another or remain distinct.1.1 Objectives of the thesisThe overall objective of this thesis is to develop models to predict ?spread-ing? of the solids (i.e. proppant) fraction during the fracturing operation. Inparticular, we focus at transport of proppant along a vertical uniform pipe11.1. Objectives of the thesis a) b)Figure 1.1: a) Schematic of hydraulic fracturing from http://en.skifergas.dk/ ;b) Particulate fluid cycle in the pipe/wellbore.(simulating transport to the fracture). We first examine typical flow regimesover a range of process parameters, in order to characterize the transport phe-nomena that are likely to be dominant. Next we develop a 1D model thatallows us to estimate: (i) dispersion of solid particles along a vertical pipe ina fully turbulent flow of a shear thinning yield stress fluid, and (ii) slip of thesolid phase relative to the mean flow. In dimensionless form, we end up with21.2. Literature overviewthe following advection-diffusion equation:????t +?F?z =??z [D????z ]. (1.1)In this equation ??(z, t) represents the average concentration of solids, z is adimensionless length, scaled with the pipe length L? and t is dimensionless time,scaled with L?/U?0, where U?0 is the mean axial velocity. The function F includesthe averaged effect of settling relative to the mean flow, and D represents thetotal diffusivity/disspersivity of the flow. As well as deriving (1.1), we giveanalytical estimates for both the advective and diffusive contributions to thespreading of solids. Finally, we implement a robust numerical scheme forsolving (1.1), which allows us to investigate various pumping strategies.1.2 Literature overviewIn general, models of fracturing hydraulics have focused on the fracture itself.Models range from simplified hydraulic descriptions to those that consider two-phase governing equations, assuming that solid (proppant) and fluid phasescan be described as two phases of incompressible continua. Here, we reviewonly models of the fracturing flows in the pipe/wellbore.The majority of fracturing flows within the wellbore are highly turbulentdue to the high flow rates and relatively small pipe diameters 1 Thus, frequentlythese are considered as fully homogeneous mixtures. In the case that the flowrates are moderate/low and the fracking fluid rheology is high, the particle-scale flows are non-inertial and appropriate models would consider the mixtureas a (non-Newtonian) viscous suspension. For higher flow rates and lowerviscosities, viscous effects are less important on the particle scale and the twophases must be considered separately.1Here, we assume that the fracturing fluid is shear thinning. Therefore, shear rateincreases with decreasing the pipe diameters. This gives smaller viscosity and consequentlyhigher Reynolds number see (1.6).31.2. Literature overviewModels for pipe flow of homogeneous Newtonian suspension abound, inparticular being developed for the mining industry, e.g. [35, 37]. These areoften however focused at critical velocities and the onset of beds in near-horizontal pipelines. For vertical pipes, the solids distribution across the pipemay be assumed to be fairly uniform, so that a 1D hydraulic approach makessense.Our focus is however on dispersion of proppant slugs along the pipe. Ifthe suspension is really treated as a single fluid mixture, then this logicallyleads to descriptions of the particle diffusivity that are similar to those for theliquid phase. There are some developments of this idea based on dimensionalanalysis coupled to data-fitting, e.g. [48]. The main point to note however isthat, as with turbulent flow of liquids, the dominant process in spreading thesolid particles axially is dispersion by the axial velocity field.Although there is considerable work on modelling of Newtonian suspen-sions, i.e. those in which the liquid carrier phase is Newtonian, there is less onnon-Newtonian suspensions. This understanding is however evolving rapidly.At the level of continuum modelling of suspension rheology, perhaps the mostcomprehensive studies that concern yield stress fluids have been carried outby Ovarlez and co-workers, [3, 5, 26, 31, 32, 47], who have used a mix of ex-perimental and modelling techniques in developing a general framework thathas been validated in some limits. We will adopt a rheological model that fitswith this general framework.For more inertial suspension flows (on a particle scale) we are still in thesituation that the dominant spreading mechanism will be through axial dis-persion, but the particle diffusivity becomes distinct from that of the liquidphase. Eskin and co-workers have developed a ?Kolmogoroff approach? toslurry pipeline flows; see [14, 15]. This approach assumes the spectral energydensity distribution of eddies, but modifies the range of the spectrum accordingto particle size/separation considerations. Essentially the small scale cut-offis increased by the presence of particles. It is considered that eddies with size41.3. Dimensional analysisbelow this cut-off transform their energy into small-scale chaotic motions ofparticles. Eskin [14] shows that the neglect of this part of the spectrum hasonly a minimal effect on the rms eddy velocity scale. In [15] the fraction ofenergy dissipated in solid-liquid interactions is estimated and is used to modifythe solids phase diffusivity.More recently, a simpler semi-empirical approach has been adopted [17],that resembles that of [48] in coupling dimensional analysis to data-fitting.This leads to a relatively simple expression for the solids diffusivity that in-cludes the effects of particle diameter ratio, ?p.1.3 Dimensional analysisWe aim to model the flow of a proppant-laden shear thinning fracturing fluid,possibly with yield stress, that is being pumped in turbulent flow downwardalong a (vertical) pipe. The flow described will depend on at least the followingset of physical parameters.? The pipe diameter, D?.? The pipe length, L?.? A representative particle diameter for the proppant (solids phase), d?p.? The liquid phase density, ??f .? The solids phase density, ??s.? Gravitational acceleration, g?.? The flow rate of the slurry, Q?, measured positive in the downwards di-rection along the pipe.? The yield stress of the particle free fluid, ??Y 0.? The consistency coefficient of the particle free fluid, K?0.51.3. Dimensional analysis? The power law index of the particle free fluid, n.? The volume fraction of solid phase, ?.The last two parameters are dimensionless. The first 9 parameters dependon a 3 independent dimensions. Throughout this thesis we will adopt thepractice of denoting all dimensional quantities by the ?? symbol. Thus, the flowis minimally described by 6 dimensionless groups, plus the volume fractionof proppant (solid phase) ?, and the power law index n. We will adopt thefollowing 6 principal dimensionless groupss = ??s??f= Density ratio, (1.2)Fr =?U?20g?D?(s? 1)= Densimetric Froude number =Imposed velocityBuoyancy velocity,(1.3)? = D?L?= Aspect ratio of the pipe, (1.4)?p =d?pD?= Particle to pipe diameter ratio, (1.5)Re = ??f U?20K?0(U?/(D?/2))n= Reynolds number =Inertial stressesViscous stresses, (1.6)B = ??Y 0K?0(U?0/(D?/2))n= Bingham number =Yield stressViscous stress in the flow.(1.7)Where U?0 = 4Q?/(?D?2) is the mean velocity. In describing different physicalphenomena and closure models, it is often convenient to use other variables,which we will do. However, these other groups may be expressed algebraicallyas functions of the above.Table 1.1 and Table 1.2 show the ranges of dimensional and non-dimensionalparameters respectively that we consider representative of the industrial appli-cation. The base set of parameters leading to the above dimensionless groups is61.3. Dimensional analysisIndustrial DataParameter Unit RangeD? (mm) 50? 120L? (m) 500? 3000d?p (mm) 0.1? 2??f (kg/m3) 1000??s (kg/m3) 2650? 3650U?0 (m/s) 2? 25??Y 0 (Pa) 0? 10K?0 (Pa.sn) 0.01? 5.0n - 0.2? 1? - 5%? 40%Table 1.1: Typical ranges of the dimensional process parameters in the indus-trial application.Dimensionless parameters Ranges 2.65? 3.65Fr 1.1? 27.8? 10?4 ? 10?3?p 0.0008? 0.04Re 3? 103 ? 107B 0? 100Table 1.2: The range of non-dimensional parameters in the industrial applica-tion.clearly not sufficient to fully describe all phenomena one is likely to encounterin a pipe flow of a turbulent suspension. Characterizing the particle distribu-tion via a single parameter is a gross simplification, (but perhaps necessary),apart from a size distribution of particles other mechanical and geometricparameters may become important, e.g. friction coefficient, maximal packingfraction, wall roughness etc...71.3. Dimensional analysisBased on the ranges of physical parameters likely to be encountered, wemay begin to characterise the flow along the pipe, into the well. This char-acterisation has been performed in [18], considering both the bulk flow andparticle scale effects. The first conclusion is that the bulk flow of most frac-turing flows when traveling down the well is a fully turbulent suspension flow.Apart from depletion layers close to the pipe walls we can expect turbulenteddies to make the leading order solids concentration uniform across the pipe.Only in more extreme cases of low flow rates and very viscous fracturing fluidswould we begin to enter a weakly turbulent or laminar regime, which wouldresult in increased dispersion relative to the mean flow. Secondly, on the parti-cle scale it was found that flows range from those where particle dynamics areessentially Stokesian (low to moderate flow rates and moderate fracking fluidrheologies) through to those where the particle dynamics may be inertial andthe interphase coupling relatively weak. The latter results from a combinationof low rheology and high shear rates.Some idea of the range of sedimentation velocities can be gained from ap-plication of classical correlations for hindered settling, e.g. Richardson andZaki [36], assuming the typical ranges for the effective viscosity given in Ta-ble 1.1. We can show that the particle settling velocity, v?p, is typically lessthan 10?3U?0. The ratio between particle settling and bulk velocities indicatesthe distance settled during transit along the pipe, compared to the length ofpipe. While this distance can be of the order of 1?3m it is important to notethat an entire slug of proppant would be settling at similar rates. Insofar asdispersion of a slug is concerned, advective dispersion is governed by the dif-ference in particle velocity as ? ranges over [0, ?in]. Taking a typical ?in ? 0.3,typically v?p could range from v?p0 to 0.25v?p0, meaning that the dispersive effectis comparable to the net settling effect. Here v?p0 represents the unhinderedsettling velocity. Therefore, it would be unwise to neglect particle settling evenin vertical pipes. Although we expect that other effects will be dominant, itis conceivable that with shorter time intervals between pulsed slugs and for81.4. Outline of the thesislong wells, settling due to advection could become relevant. Note also that ad-vective spreading is linear in transit time (t?transit = L?/U?0), whereas diffusivespreading should increase like t?1/2transit.1.4 Outline of the thesisIn chapter 2, we derive a model for particle transport, reduce the model viaphysical scaling arguments and finally use the method of multiple-scale to de-rive a 1D advection-diffusion model for transport and dispersion of the meansolid particle concentration along the wellbore. The advection part of themodel includes the transport of the solid phase with the mean flow velocityand the mean particle settling velocity. The diffusion part of the model in-cludes both the averaged turbulent particle diffusivity and the effect of Taylordispersion. In chapter 3, we develop a robust numerical algorithm to solvethe nonlinear 1D advection/diffusion model. This algorithm is used to studythe spreading of the slug of solid along the well for few industrial examples.The thesis concludes in chapter 4 with a summary and number of generalrecommendations.9Chapter 2Modeling approachThe overall objective of this chapter is to derive a 1D advection diffusionequation (1.1) to estimate dispersion of solid particles along a vertical pipein a fully turbulent flow of a shear thinning yield stress fluid. The advectionpart of the model includes the transport of the solid phase with the meanflow velocity and the mean particle settling velocity. The diffusion part of themodel includes both the averaged turbulent particle diffusivity and the effectof Taylor dispersion. We also provide a simple analysis of (1.1), suitable forspreadsheet-type field design purposes, in which we estimate ?mixing lengths?due to to both settling and dispersion.2.1 FormulationWe start with the governing equations for two incompressible continuous phases.Then, we derive the governing equations for the mixture and explain the trans-port equation for the dispersed phase. We explain the constitutive law for themixture and we obtain the averaged velocity profile for a fully developed tur-bulent flow of a shear thinning yield stress fluid in a pipe. Finally, we give theclosures for the solid-fluid interaction force and relative velocity.2.1.1 Mass conservation and particle diffusivityWe assume that solid and fluid phases can be described as two phases ofincompressible continua, which implies some form volume-averaging over ascale larger than the particle scale in interpreting the flow variables. The mass102.1. Formulationconservation equations for each phase are:??s???t?+ ?? ? (??s?u?p) = 0, (2.1)??f?(1? ?)?t?+ ?? ? (??f (1? ?)u?f) = 0. (2.2)Where ? is the local volume fraction of solid. The phase-averaged solid andfluid velocities are denoted by u?p and u?f respectively. On dividing through(2.1) by ??s and (2.2) by ??f , then summing:?? ? u? = 0 : u? = ?u?p + (1? ?)u?f . (2.3)Here u? is the volume-averaged velocity for the mixture of solid and fluid phases.Equation (2.1) can be rewritten, in terms of u? and the relative velocity betweensolid and liquid phases:???t?+ ?? ? [?u?+ ?(1? ?)u?r] = 0, (2.4)where u?r = u?p ? u?f .In practice, as we consider turbulent flows, we are more interested in theaverage behaviour of the flow variables, which may evolve temporally andspatially. Therefore, we consider ensemble-averaged quantities, denoted withan over-bar, and fluctuating components:? = ??+ ??, u?p = ??up + u??p, , u?f = ??uf + u??f , etc.Ensemble averaging of (2.4) leads to:????t?+ ?? ? [????u+ ??(1? ??)??ur] = ??c, (2.5)where the term ??c results from correlation of the fluctuating components of112.1. Formulationthe velocity and solids fraction fields. More formally, we can identify ??c:??c = ??? ? [??u?? + (1? 2??)??u??r ? (??)2 ??ur] ? ??? ? [??u??], (2.6)on assuming that the relative velocity terms as relatively small. Typically now,the quantity ??c is modeled as a diffusive flux via Fick?s law, in order to closethe equations, i.e.??u??p = ?D?d????, (2.7)where D?d is the particle diffusivity coefficient. An expression for D?d has to bedetermined either from experiment or through some auxiliary analysis.Finally, we ensemble-average (2.3) to give:?? ? ??u = 0, (2.8)and the solids phase mass conservation approximation:????t?+ ?? ? [????u+ ??(1? ??)??ur] = ?? ? [D?d????], (2.9)We note that (2.9) contains the main components necessary to model disper-sion along the pipe. It is necessary to provide closure expressions for the meanaxial velocity, the relative velocity and the particle diffusivity.2.1.2 Momentum balanceThe linear momentum balances for each phase can be written as:??s?(?u?p?t?+ (u?p ? ??)u?p)= ?? ? ??p + ??s?g?k + m?, (2.10)??f (1? ?)(?u?f?t?+ (u?f ? ??)u?f)= ???p?f + ?? ? ??f + ??f(1? ?)g?k ? m?.(2.11)122.1. FormulationWhere ??p denotes the particle stress tensor, which may be further decomposedinto shear and normal stress components. The fluid phase shear stress tensoris denoted ??f and m? denotes the solid-liquid interaction force, per unit mass.Normally however, we are concerned with the bulk mixture momentumequation, written as the sum of the (2.10) and (2.11). We also assume thatthe flow varies slowly in the axial direction, due to variations in mean solidsconcentration. For a fully developed turbulent flow, the averaged componentsof the inertial components on the left hand side of (2.10) and (2.11) vanish,leaving only Reynolds stress terms on the the right-hand side. The solid-liquidinteraction terms cancel out. We write this simplified mixture momentumbalance as:0 = ???p?+ ?? ? ?? + ??g?k, (2.12)where we decompose the shear stress tensor ?? into a viscous term and a tur-bulent term, i.e.?? = ?? v + ?? t. (2.13)2.1.3 Viscous stress closureThe approach adopted for modelling the viscous stress tensor is outlined in[18]. The fracking fluid is assumed to be of Herschel-Bulkley type. We fol-low the developments of Ovarlez and co-workers, [3, 5, 26, 31, 32, 47], whohave used a mix of experimental and modelling techniques in developing ageneral framework that has been validated in some limits for non-Newtoniansuspension flows. The general framework is as follows.Shear flows of suspensions are characterised by a bulk suspension viscosity??, which relates to the viscosity of the mixture of solid and liquid phasesand is used to model the viscous component of the shear stress tensor in aconventional way, i.e.?? v = ?? ???(??u), (2.14)where ???(??u) is the rate of strain tensor. The suspension viscosity is first de-132.1. Formulationcomposed as follows?? = ??f?r(??), (2.15)where ??f is referred to as the liquid phase viscosity and ?r(??) is the dimension-less relative viscosity. We will assume dependency on the ensemble averagedsolid fraction ?? rather than on ?, as we wish to work with the averaged quan-tities. This decomposition of the suspension viscosity is quite classical and therelative viscosity is modeled by a closure law, e.g. most simply the Einstein-Roscoe law. Since here we deal with significant ??, we suggest adopting theKrieger-Dougherty law:?r(??) =[1? ???m]?2.5?m, (2.16)or one of its close variants, e.g. [12, 21, 27, 34]. Apart form conceptual simplic-ity, there exist generalisations to particles of different shapes e.g. rods/fibres,see [49]. For the results presented in this thesis we will assume a maximalpacking fraction of ?m = 0.57.The relative viscosity can be further decomposed to model the solid phasestress tensor, but this is not needed here. For the liquid phase viscosity ??f , twoeffects must be considered. Firstly, the fluids are shear-thinning and secondlythe particles modify the viscosity. In the absence of particles the effectiveviscosity is given by a constitutive law that depends on the rate of strain,i.e. the pure fracking fluid has effective viscosity:??f,0(???loc) = ??0 ???n?1loc +??Y 0???loc, (2.17)where ??0, ??Y 0 and n are the consistency, the yield stress and the power lawindex of the pure fracking fluid, respectively. Here ???loc denotes the (local)strain rate of the fluid.A significant contribution of Ovarlez and co-workers is in recognising thatalthough the solids phase increases the viscosity of the suspension via ?r(??), the142.1. Formulationpresence of particles also reduces the viscosity of the inter-particle fluid, whenit is shear-thinning. We denote the bulk suspension strain rate by ???, which iscomputed as the second invariant of the tensor ???(??u). The main point is that??? is effectively amplified by the presence of particles, since velocity gradientsare concentrated only within the liquid phase. Ovarlez and co-workers havedeveloped the following law to describe this amplification:???loc =[?r(??)1? ??]1/2???, (2.18)which is verified reasonably well by experimental results [26, 31] and also agreeswith theoretical considerations, [3].This leads to the following definition of the liquid phase viscosity ??f :??f(???, ??) = ??f,0(???loc(??)) = ??0[?r(??)1? ??](n?1)/2???n?1 + ??Y 0???[?r(??)1? ??]?1/2. (2.19)This can also be interpreted as assuming that the liquid within the suspensionobeys a Herschel-Bulkley type law, but with ??-dependent consistency and yieldstress defined by:?? = ??0[?r(??)1? ??](n?1)/2??Y =??Y 0[?r(??)1???]1/2 . (2.20)The framework given by (2.15)-(2.21) is quite simple to work with and modify.Alternatively, we can include directly the bulk viscousification effect of therelative viscosity to get:??s(??) = ??0Y (??) ??Y s(??) = ??Y 0X(??) (2.21)X(??) =?(1? ??)?r(??) Y (??) = ?r(??)[?r(??)1? ??](n?1)/2152.1. Formulationa) 0 0.1 0.2 0.3 0.4 0.501234X??b) 0 0.1 0.2 0.3 0.4 0.505101520Y??n=1n=0.2Figure 2.1: a) X vs ??; b) Y vs ?? for n = 0.2, 0.4, 0.6, 0.8, 1.0.Therefore, we can write the effective suspension viscosity as follows??(???, ??) = ??s[?r(??)1? ??](n?1)/2???n?1 + ??Y s???[?r(??)1? ??]?1/2. (2.22)Figure 2.1 plots the functions X(??) and Y (??) that describe the amplifi-cation of the yield stress and consistency within the suspension. In all caseswe see a viscousification effect as ?? increases. Note however, that this is as aresult of the relative viscosity: on neglecting the final multiplication by ?r(??),the functions in (2.21) are decreasing with ?? ,i.e., local viscosity on particlescale (2.17) is decreasing with ??, due to shear rate amplification by particles.2.1.4 Turbulent velocity closureWe consider fully turbulent flow along a vertical pipe. The ensemble-averagedsolids fraction ?? is assumed uniform across the pipe. The following system of162.1. Formulationequations governs the flow of the mixture:0 = ??p??z? +1r???r? (r???rz) + ??g?, (2.23)0 = ??p??r? +1r???r? (r???rr), (2.24)0 =1r???r? (r???u) + ???w?z? , (2.25)where ??rz and ??rr can be decomposed to viscous and turbulent parts. Theviscous component has been discussed above. The turbulent component arisesfrom Reynolds stress terms, approximated as follows:??rz = ?? vrz + ?? trz = ?? vrz ? ??(u??w??), (2.26)??rr = ?? vrr + ?? trr = ?? vrr ? ??(u??u??). (2.27)We assume that the Reynolds stress (?? trr) is independent of z and themodified pressure gradient?P??z? =?p??z? ? ??g?,is constant across the pipe cross-section. We integrate equation (2.23) fromwall of the pipe y? = 0 outwards:?? trz + ?? vrz = ??w +y?2?P??z? : ??w = ?R?2?P??z? . (2.28)Note that ??w is the wall shear stress and y? is R? ? r?. This is used to define aturbulent velocity scale W? ?:?? vw = ??(W? ?)2. (2.29)172.1. FormulationWe substitute (2.29) into (2.28) to give?? trz + ?? vrz = ??(W? ?)2(1? y?R?). (2.30)The viscous stress is mainly large in the vicinity of the wall, where ?? trz is small.According to [38], the proper scaling factor for ?? vrz is:?? vrz = ??w? ??w?y? : ??w =??w[(??w ? ??Y s)/??s]1/n. (2.31)In the core region of the pipe where Reynolds stress is dominant, we definey = y?/R? and scale the Reynolds stress: ? trz = ?? trz/[??(W? ?)2], to give:? trz +1Re???y(??wW? ?)= (1? y). (2.32)The Reynolds number Re? is defined by:Re? = ??W??R???w.We typically expect that Re? ? 1 when fully turbulent, meaning that theviscous stresses are largely irrelevant in the core.On the other hand, the stress at the pipe wall is purely viscous, so thatthe scaling for (2.32) can not be valid near the wall in the limit Re? ? ?.From (2.32) we conclude that the near-wall behaviour can be accounted for byabsorbing Re? in the scaling for y?, i.e. y+ = yRe? = y?W? ???/??w. The resultingequation is? trz +??y+(??wW? ?)= (1? y+Re? ). (2.33)When the production and dissipation of the mixture turbulent energy are182.1. Formulationequal (i.e. in a state of local equilibrium, see [43]), we may write? < u??w?? > (???w?y? ) ?W ?3R?as Re? ?? : ?< u??w?? >W ?2 = (1? y)? ???w?y? =W ?R??F (y)?y . (2.34)Equation (2.33) implies that w?/W ? = G(y+) and consequently ? ??w/?y? =(W ?Re?/R?)(?G/?y+). In the limit that y+ ?? and y ? 0, we may writeW ?R?F (y)?y =W ?Re?R?G?y+ ?y+n?G?y+ =yn?F?y =1?n. (2.35)Finally, we have??wW? ?=1?1nln y+n + ?2n . (2.36)Note that the log-law is valid in the matching region where y+ ? ? andy ? 0. However, this is an engineering approximation to consider the log-law velocity profile for the whole cross section of the pipe. Following [38],who study weakly turbulent power law and Herschel-Bulkley fluids using DNStechniques, the suggested values for the constants in (2.36) are ?1 = 0.4 and?2 = 5.5.We assume that the dimensionless flow rate across the pipe is ? in orderto calculate W ?. We define the friction factor in the usual way:ff = 2(W? ?/U?0)2,where U?0 is the mean velocity of the mixture. Then we scale w? = ??w/U?0 and192.1. Formulation0.5 1 1.5 2 2.5x 105510152025U?0(m/s)Rea)0.2 0.4 0.6 0.8 1 1.2510152025U?0(m/s)Bb)5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 0.002 0.004 0.006 0.008 0.01c)Figure 2.2: a) & b) Dependency of bulk Re and B on mean flow; c) colormapof friction factor in (U?0, ??) plane (see equation 2.38). For all figures, we haven = 0.5, K?0 = 0.1Pa.sn, ??Y 0 = 1Pa, s = 2.65, D? = 62mm.write (2.36) in the following dimensionless form.w? =?0.5ff(?2n +1?1ln[(1? r)Rew(0.5ff)1/2]), (2.37)where r = r?/R? andRew =??R?U?0??w.202.1. Formulation1 2 3 4 5 6 7x 1050.20.61nRea)0.5 1 1.5 2 2.5 30.20.61nBb)0.2 0.4 0.6 0.8 100.10.20.30.40.5n?? 0 0.01 0.02 0.03 0.04 0.05c)Figure 2.3: a) & b) Dependency of bulk Re and B on power law index n; c)colormap of friction factor in (n, ??) plane (see equation 2.38). For all figures,we have U?0 = 15m/s, K?0 = 0.1Pa.sn, ??Y 0 = 1Pa, s = 2.65, D? = 62mm.In order for the flow to be fully turbulent, it is expected that Rew ? 1. Tocompute the friction factor, we integrate across the pipe cross-section, usingthe constraint that the flow rate is ?, or the mean velocity is unity:1 =? 102rw?(r) dr.212.1. FormulationThis gives us, at leading order:1 =?0.5ff(?2n +1?1ln[Rew(0.5ff)1/2]? 1.5?1), (2.38)We obtain the dependency of Rew on the dimensionless parameters of the bulkflow as follows.Rew =Re1/n0.5ff[(0.5ffY (??))(s?+ 1? ?)?( BReX(??)Y (??))]1/n. (2.39)Equation (2.39) is obtained by expressing mean flow velocity and stress in(2.38) in terms of bulk Re and B. Where, B is ??Y 0/??0(U?0/R?)n and Re is??f U?20 /??0(U?0/R?)n. Two additional important effects that one might try to addare: (i) some form of depletion of the particle phase near the wall; (ii) wallroughness effects. Both these effects are however most likely to be significantwith respect to dispersion only for weakly turbulent flows, when the wall layersbecome thicker.Figures 2.2a & b show the the dependency of bulk Re and B on the velocityof mean flow U?0 respectively. Re increases with the velocity of mean flowbecause both inertia and shear rate increase in the system. The fluid is shearthinning the latter decreases the viscosity. The total stress increases with thevelocity of mean flow and consequently B decreases. Figure 2.2c shows thecolormap of ff in the (?, U?0) plane. We see that the friction factor increaseswith the solids volume fraction as a result of viscousification. Moreover, thefriction factor decreases with velocity of mean flow. Although total stressincreases with U?0, the ratio of W? ?/U?0 decreases with U?0.Figure 2.3 shows similar results, but with variations in n rather than U?0.Figures 2.3a & b show that the bulk Re and B decrease with n. The reasonis that the viscous stress increases with n. Figure 2.3c shows that the frictionfactor increases with both the solids volume fraction and n as a result ofviscousification.222.1. Formulation2.1.5 Relative velocity closureIn this section, we develop the closure for m? in terms of the relative velocitybetween solid and fluid phases, ??ur, and ??. Principally this involves hydrody-namic effects of viscous drag:m? ? F?D6???d?3p. (2.40)whereF?D = ???f2|??ur|??ur(?d?2p/4)CD. (2.41)The interaction term is also estimated from combining (2.10) and (2.11):m? = ??(1??)[??s???f ]g?k+?(1??)[??sDpu?pDt?? ??fDf u?fDt?]?(1??)?????p+?[???p?f+?????f ].(2.42)It is usual to assume that the first term is dominant, leading to:|??ur|??ur =43CD(1? ?)(s? 1)g?d?p (2.43)Despite numerous studies on the motion of a particle in Newtonian fluids,there are only a limited number of studies on the motion of an object in theyield stress fluids. As well as non-Newtonian effects, we need to incorporatethe effect of the particles on the settling, i.e. via hindering. The approach thatwe adopt follows that of [32]. We compute the single particle settling speedfrom (2.43). Similar to a Newtonian fluid we adopt a drag coefficient closuresuch as:CD =24Rep(1 + 0.1Re0.75p ); (2.44)see e.g. [20]. Here Rep = ??f |??ur|d?p/??f . Recall that ??f is the liquid phaseviscosity around the particles; see (2.19). Combining (2.44) with (2.43) leadsto a unique determination of the relative velocity for a single particle, say ??ur,0.232.1. FormulationFor example, assuming NRe ? 1 we have the Stokes solution:??ur,0 =118(1? ?)(s? 1)??f g?d?2p??f. (2.45)The single particle settling velocity is now multiplied by a hindering functionh(??), to give the relative velocity of the suspension, i.e.??ur = h(??)??ur,0. (2.46)Following [32], we useh(?) = 1? ??r(?). (2.47)A similar approach has been adopted by [41].It remains to specify the suspension shear rate within the turbulent flow,which is needed to evaluate the effective viscosity ??f . Close to the walls, thevelocity changes from U?0 to zero over approximately the length-scale of theviscous sublayer, i.e. an effective shear rate would be something like:??? = |??ur|d?p+U?0y?sub(2.48)However, at high Reynolds numbers the viscous the wall region has only aminor effect on dispersion, compared with the turbulent core. In the core,there are both changes in mean flow (with shear rate roughly ? W? ?D?), and inthe fluctuating component of the velocity. A very simple model is to assumethe velocity fluctuations of size ? W? ? occur on a scale of the largest eddies,assumed to have size roughly 0.1D?. This leads to:??? ? |??ur|d?p+W? ?0.1D?. (2.49)Note that this model suggests that the strain rates due to the velocity fluctu-242.2. Particle dispersionations are probably the dominant component in the turbulent core (comparedto gradients in the mean velocity), which has the consequence that the relativevelocity will vary with r? only due to variations in ??.Figure 2.4 shows colourmaps of various flow parameters in the (U?0, ??) plane.We note that the total shear rate increases with both U?0 and ??. The formeris due to increase in bulk shear rate and the latter is as a result of shear rateamplification by particles. We calculated ??? for a range of flow parameters andobserved that the main contribution to ??? comes from background shear rate,i.e. the second term in (2.49). Since the suspension is shear thinning, ??(???) is adecreasing function of ???. On the particle scale, we see that the relative velocity,??ur, decreases with viscousification (i.e. the increase in ??) and increases withliquefaction (i.e. the increase in U?0 and bulk shear rate). Therefore, we seethat Rep decreases with ?? and increases with U?0.Figure 2.5 shows colourmaps of various flow parameters in the (d?p, ??) plane.Note that the friction factor (2.38) and hence W? ? do not depend on the sizeof particles, d?p. Since, ??? mainly originates from the background shear rate,i.e. W? ?/0.1D?, we can conclude that the total shear rate, ???, and suspensionviscosity, ??(???), will have only a small variation with d?p. The settling velocityof particles does however increase with the particle size due to the volumetricincrease in gravity force. Therefore, ??ur and consequently Rep are increasingfunctions of the particle size.2.2 Particle dispersionIn most pipeline transport processes, spreading of a concentration (e.g. solidsfraction), is governed by a number of contributions, from both diffusive andadvective processes. In the case of fully turbulent flows in moderately longpipes, it is possible to estimate the net effect of the advective component inspreading concentrations axially, in what has been termed Taylor dispersion[42]. The net effect of axial spreading by advection results in diffusion of the252.2. Particle dispersion5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 50 100 150a)???5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 0.1 0.2 0.3 0.4 0.5b)??(???)5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 0.1 0.2 0.3 0.4 0.5 0.6c)Rep5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 0.001 0.002d)??urU?0Figure 2.4: Colourmaps of various flow parameters plotted in the (U?0, ??) plane:a) total shear rate, ???(1/s), as calculated from equation (2.49); b) effectiveviscosity of the suspension, ??(???)(Pa.s), as calculated from (2.22); c) particle Renumber, Rep; d) dimensionless slip velocity ??ur/U?0, as calculated from (2.43).For all figures, we have n = 0.5, K?0 = 0.1Pa.sn, ??Y 0 = 1Pa, s = 2.65, ?p = 0.01,D? = 62mm.262.2. Particle dispersion4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 0 10 20 30 40 50a)??urd?p4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 90 95 100 105 110 115 120b)W??0.1D?4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 0 1 2 3 4 5 6c)Rep4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 0 0.002 0.004d)??urU?0Figure 2.5: Colourmaps of various flow parameters plotted in the (dp, ??) plane:a) particle shear rate ??urd?p ; b) background shear rateW? ?0.1D? ; c) particle Reynoldsnumber, Rep; d) dimensionless slip velocity ??ur/U?0, as calculated from (2.43).For all figures, we have U?0 = 15m/s, n = 0.5, K?0 = 0.1Pa.sn, ??Y 0 = 1Pa,s = 2.65, D? = 62mm.272.2. Particle dispersionmean concentration, relative to a constant advective speed. Importantly, theTaylor dispersion effect dominates diffusivity effects. In this section, belowwe develop the general expression for the Taylor dispersion effect, applied tothe solids fraction ??. Later we compute the net diffusive effects for 3 differentparticle diffusivity models, and give estimates for the length of pipe over whicha slug of proppant would disperse in travelling down the pipe.2.2.1 Dispersion modelFollowing the analysis of ?2.1.4, the mean axial velocity in the turbulent coreis given by:??w = W? ?(1?1nln y+n + ?2n), (2.50)which is coupled to the mean solids fraction ?? through the rheological closuresentering the wall-coordinate scaling for y+. The mean velocity and solidsfraction satisfy:0 =1r???r? (r???u) + ???w?z? , (2.51)????t?+ ??u????r? +??w????z? +??????????z? =1r???r? (r?D?d????r? ) +??z? (D?d????z? ), (2.52)where ?? = ??(1? ??) ??wr.We assume that D?d = D?d(r?, ??) and write:D?d(r?, ??) = D?d,aveDd(r?, ??) : D?d,ave =2?max? ?max0? R?0D?d(r?, ??)r? dr?d??,with ?max some upper limit for ??. We use the averaged particle diffusivity to282.2. Particle dispersiondefine a Pe?clet number Pe via:Pe = R?U?0D?d,ave. (2.53)Note that typically the Pe?clet number is larger than unity (approximatelyscaling with a negative fractional power of the friction factor, depending onthe model used - see below). The Pe?clet number can be interpreted as thenumber of radii that the mixture needs to advect at mean speed in orderfor the particle diffusivity to diffuse the solids concentration across the pipe.Although greater than unity, it is also common that ?Pe ? 1, since pipeaspect ratios are relatively long.Equations (2.51) , (2.52) and the velocity profile are made dimensionlesswith the following scaling:z = z?L?, r = r?R?, t = t? U?0L?, u? = 2??u?U?0, w? =??wU?0, w?r =??wrU?0, ? = ??U?0,?0.5ff =W? ?U?0.On substituting the non-dimensional variables into the governing equations,we arrive at:w? =?0.5ff(?2n +1?1ln[(1? r)Rew(0.5ff)1/2]), (2.54)0 =1r??r (ru?) +?w??z , (2.55)?Pe2[????t + u?????r + w?????z +?????????z]=1r??r (rDd????r ) +?24??z(Dd????z)(2.56)Recall that ? = D?/L? = 2R?/L?.Note that by construction, Dd(r, ??) ? 1 and the average of the mean axialvelocity w? across the pipe is 1. We now define the averaged and fluctuating292.2. Particle dispersionsettling rates:w?r,ave =? 102r??(??)??? dr,??wr =??(??)??? ? w?r,ave,and move to a moving frame of reference, by writing:??w = w? ? 1, ? = z ? (1 + w?r,ave)t.Further, we introduce ? = ?Pe/2 and infer from (2.56) that as well as theadvective timescale (captured in t), the mean concentration will also respondon the slow timescale: T = ?t. We do the latter by a multiple timescalesapproach, assuming that ?? = ??(t, T, r, ?).Equations (2.55) and (2.56) become, in the moving frame0 =1r??r (ru?) +? ??w?? , (2.57)?[????t + ?????T + u?????r +??w????? +??wr?????]=1r??r (rDd????r )+?2Pe2???(Dd?????). (2.58)These are coupled with the following boundary/symmetry conditions????r = 0 at r = 0, 1, (2.59)u? = 0 at r = 0, 1. (2.60)We now exploit the fact that ? ? 1 to derive a leading order approximationto the behaviour of ??, using the following asymptotic expansion??(t, T, r, ?) = ??0(t, T, r, ?) + ???1(t, T, r, ?) + ?2??2(t, T, r, ?) + ? ? ? ,u?(t, T, r, ?) = u?0(t, T, r, ?) + ?u?1(t, T, r, ?) + ?2u?2(t, T, r, ?) + ? ? ? ,??w(t, T, r, ?) = ??w0(t, T, r, ?) + ? ??w1(t, T, r, ?) + ?2 ??w2(t, T, r, ?) + ? ? ? .302.2. Particle dispersionZero-th order problem:At leading order we have:0 =1r??r (ru?0) +? ??w0?? , (2.61)0 =1r??r (rDd,0???0?r ). (2.62)???0?r = 0 at r = 0, 1,where the additional subscripts mean that the expressions are expanded withrespect to the leading order variables. From (2.62) we see that ?0 has no r?-dependency. One consequence of this is that, since the relative velocity variesradially only via ??, we have that ? = ??(1? ??)w?r is also independent of r andhence:??wr,0 = 0.First order problem: At first order we have:???0?t +??w0???0?? =1r??r(rDd,0???1?r), (2.63)???1?r = 0 at r = 0, 1.with all other terms vanishing since ?0 has no r?-dependency. We then multiplyby r and integrate (2.63) with respect to r over [0, 1]. The right-hand sidevanishes due to the boundary conditions on ?1. On the LHS, we find only:???0?t = 0, as ??w0 has zero average. It follows that at leading order ??0 = ??0(T, ?).Returning to (2.63), we now find:??1(t, T, r, ?) = ??1(t, T, 0, ?) +???0?? (T, ?)?(T, r, ?), (2.64)312.2. Particle dispersionwhere?(T, r, ?) =? r01sDd0(? s0y ??w0(T, y, ?) dy)ds. (2.65)Second order problem: At O(?2), we find that ??2 satisfies:???0?T +???1?t + u?0???1?r +??w0???1?? +??w1???0?? +??wr,1???0?? =1r??r(rDd,0???2?r + rDd,1???1?r)+1Pe2???(Dd0???0??), (2.66)???2?r = 0 at r = 0, 1.We multiply by 2r and integrate (2.66) with respect to r over [0, 1]. We use(2.61) to combine the 3rd and 4th terms on the left hand side. After theintegration, only the following terms are non-zero:???0?T +???1?t (t, T, 0, ?) +???? 102r??1 ??w0 dr =1Pe2?2??0??2 .Now we substitute for ??1 and integrate by parts the third term:???0?T +???1?t (t, T, 0, ?) =???([DT +1Pe2] ???0??).where DT is the Taylor dispersion coefficient, given by:DT = 2? 101sDd,0[? s0y ??w0(T, y, ?) dy]2ds. (2.67)Note that the only term that depends on t is ???1?t (t, T, 0, ?). It is not hardto deduce that ??1(t, T, 0, ?) = A(T, ?)t + B(T, ?), so that the only boundedsolution for ??1 must have A(T, ?) = 0.Synopsis: the axial diffusion equation322.2. Particle dispersionFinally, this leads to the following diffusion equation for ??0(T, ?):???0?T =???([DT +1Pe2] ???0??), (2.68)which can be written in the stationary frame of reference as:???0?T +[1 + w?r,ave]????0?z =??z([DT +1Pe2] ???0?z), (2.69)or in conservative form and in terms of the fast time variable:???0?t +??zF (??0) =??z(?[DT +1Pe2] ???0?z): F (??0) = ??0[1 + (1? ??0)w?r].(2.70)To conclude, we rescaled the problem into dimensional form:???0?t?+ [U?0 + ??wr,ave]???0?z? =??z?(D?d,ave[1 + Pe2DT] ???0?z?), (2.71)We observe that the leading order solids fraction diffuses relative to themean motion of pumping plus the averaged settling. The diffusive spreadingis amplified by Taylor dispersion, resulting in an effective axial diffusivity of(1 + Pe2DT )D?d,ave, where D?d,ave is the averaged particle dispersivity. TheTaylor dispersion part of the axial diffusivity can be written in dimensionalform as:D?T = Pe2DT D?d,ave =R?2U?20D?d,ave2? 101sDd,0[? s0y ??w0(T, y, ?) dy]2ds= 2? R?01s?D?d,0[1R?? s?0y?( ??w0 ? U?0) dy?]2ds?. (2.72)332.2. Particle dispersion2.2.2 Particle diffusivity modelsTo get an idea of the size of the axial diffusivity, we evaluate the Taylordispersion coefficient for 3 sensible particle diffusivity models. In order tomake comparisons between models, and because the dispersive part of theaxial diffusivity is dominant, we adopt the same friction factor and velocityprofile for each model.The velocity profile is given dimensionally by (2.50) and dimensionlesslyby (2.54). Both friction factor and velocity profile are derived in ?2.1.4. Wenote that ??w0 is given by??w0 =?0.5ff?1(ln(1? r) + 32). (2.73)and in the definition of DT , we can evaluate:? r0y ??w0(T, y, ?) dy = ??0.5ff2?1(1? r) (r + (1 + r) ln(1? r)) . (2.74)Homogeneous slurry/Reynolds analogy approachThis approach treats the slurry as a fluid-solid mixture, modeled locally withan effective viscosity, and then applies Reynolds analogy. The mixture mo-mentum balance along the pipe (2.23), is effectively:r?2?P??z? = ??rz = ??e? ??w?r? ? ??e = ?r?R???(W? ?)2[? ??w?r?]?1, (2.75)where ??e represents the effective viscosity. We apply Reynolds analogy toapproximate the diffusivity coefficient of the particles, D?d, which is assumedequivalent to that of the liquid in the mixture (i.e. we equate the diffusivetransport of mass and momentum by the turbulent eddies). This leads to:D?d =??e?? ? ?r?R?(W? ?)2[? ??w?r?]?1= R?W? ??1r(1? r) (2.76)342.2. Particle dispersionTherefore, we find:D?d,ave =R?W? ??16= R?U?0?1?0.5ff6(2.77)Dd,0 = 6r(1? r). (2.78)The Pe?clet number is evaluated as: Pe = 6/(?1?0.5ff) and we compute:DT (??0) =ff24?21? 10(1? r)2[r + (1 + r) ln(1? r)]2r[r(1? r)] dr ? 0.0851ff . (2.79)Finally, the dimensional total axial diffusivity is:(1 + Pe2DT )D?d,ave ? [1 + 19.155]D?d,ave ? 0.95R?U?0?ff (2.80)Note that approximately 95% of the total axial diffusivity is due to Taylordispersion.WaltonIn [48] Walton uses dimensional analysis arguments to arrive at an assumedfunctional form of the particle diffusivity closure relation:D?d = U?0R???p0 : ??p0 ? ?a1Rea2 . (2.81)The coefficients in the above expression are fitted to existing experimentaldata collected by [30]. Note that there is no radial dependence of D?d.For simplicity we take the same approach. However as discussed above,we substitute our friction factor closure from equation (2.38), in order to in-clude both non-Newtonian and solid volume fraction effects. This leads to thefollowing expression for the particle diffusivityD?d = D?d,ave = 2U?0R??0.5ff0.014?0.96p Re0.33?0.5? 0.079Re?0.25. (2.82)352.2. Particle dispersionIt follows thatDd,0 = 1, (2.83)Pe = 12?0.5ff?0.5? 0.079Re?0.250.014?0.96p Re0.33, (2.84)DT (?0) =ff4?21? 10[(1? r) (r + (1 + r) ln(1? r))]2r dr ? 0.1703ff(2.85)Finally, the dimensional axial diffusivity is:(1 + Pe2DT )D?d,ave ? [1 + 17.16Re?0.91??1.92p ]D?d,ave, (2.86)which shows a strong dependency on particle diameter ratio and bulk Reynoldsnumber.EskinEskin?s approach [17] is in some ways similar to Walton [48] , although therationale is different. In [17] the following model for D?d,ave is derived.D?d,ave = 2U?0R?? 0.292?0.32p f2/3f(1 +Rep18U?0??wr,0s?pCD0CD)(2.87)The second term above is included to account for potential particle relaxationeffects in flows that are inertial on the particle scale. Again no variation with362.2. Particle dispersionr is considered and hence we find the following expressions:Dd,0 = 1, (2.88)Pe = 10.584?0.32p f2/3f11 + Rep18U?0??wr,0 s?pCD0CD, (2.89)DT (?0) =ff4?21? 10[(1? r) (r + (1 + r) ln(1? r))]2r dr ? 0.1703ff(2.90)The dimensional axial diffusivity is:(1 + Pe2DT )D?d,ave ????1 +4.58??0.64p f?1/3f[1 + Rep18U?0??wr,0s?pCD0CD]2???D?d,ave, (2.91)Comparison of the particle diffusivity modelsFigure 2.6 shows variations in Pe for the three particle diffusivity modelsoutlined above. We plot Pe against both (U?0, ??) and (d?p, ??). We note thatthe range of Pe number is similar for both the Homogeneous slurry/Reynoldsanalogy approach and Eskin?s approach [17]. However, the homogeneous slurryapproach has no sensitivity to particle size and variations with U?0 at lowerconcentrations ?? differ from Eskin?s model. Pe is a decreasing function ofparticle size in both Walton?s approach [48] and Eskin?s approach [17]. Whichimplies that particle diffusivity increases with particle size. The Homogeneousslurry/Reynolds analogy approach does not depend on particle size.Figure 2.7 shows DTPe2 for Walton?s approach [48] and Eskin?s approach[17]. It can be seen that Taylor dispersion is the main mechanism for axialdiffusion of the mean solids fraction. Eskin?s approach [17] gives similar rangeforDTPe2 as does the Homogeneous slurry/Reynolds analogy approach, whichwe have omitted for brevity. In general we can see that for large parameterranges over 95% of the axial dispersion is due to Taylor dispersion, (see (2.91)).372.2. Particle dispersion5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 300 400 500 600a)Pe4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 400 450 500 550b)Pe5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 60 120 180 240c)Pe4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 50 100 150 200d)Pe5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 200 300 400 500e)Pe4 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 100 200 300 400 500 600 700f)PeFigure 2.6: Variations in Pe for the different particle diffusivity models: a)& b) Homogeneous slurry/Reynolds analogy approach; c) & d) Walton?s ap-proach [48]; e) & f) Eskin?s approach [17]. For all figures, we have n = 0.5,K?0 = 0.1Pa.sn, ??Y 0 = 1Pa, s = 2.65, D? = 62mm. For a), c) and e) ?p = 0.01;for b), d) and f) U?0 = 15m/s.382.2. Particle dispersion5 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 2 4 6 8 10a)DT Pe25 10 15 20 2500.10.20.30.40.5U?0 (m/s)?? 0 20 40 60 80b)DT Pe24 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 0 2 4 6 8 10c)DT Pe24 6 8 10 12 14x 10?400.10.20.30.40.5d?p (m)?? 0 20 40 60 80 100 120d)DT Pe2Figure 2.7: Variations in DTPe2 for the different particle diffusivity models;a) & c) Walton?s approach [48]; b) & d) Eskin?s approach [17]. For all figures,we have n = 0.5, K?0 = 0.1Pa.sn, ??Y 0 = 1Pa, s = 2.65, D? = 62mm. For a) andb) ?p = 0.01; for c) and d) U?0 = 15m/s.392.2. Particle dispersionWe need to be cautious about our interpretations of predictions using Wal-ton?s or Eskin?s models, which are included mainly for comparison. We seethat the majority of the solids transport is governed by axial dispersion, whichis related to the the velocity profile. Here we have used the same velocity pro-file for all models. However implicitly, the semi-empirical relations that definethe turbulent particle dispersivity in terms of a friction factor may be basedon different interpretations of the slurry velocity field.2.2.3 Slug spreading estimatesIf we observe the general form of (2.71) and consider the behaviour of a slug ofproppant pumped down the pipe at mean speed U?0, we see that there are twomain mechanisms that lead to spreading of the slug: (i) slip of the solids phasedue to settling, ahead of the mean flow; (ii) diffusion of the slug, primarilythrough Taylor dispersion. Denoting the length of well by L? the transit timeto reach the bottom of the well is L?/U?0.(i) The effect of settling is captured by the term ??wr,ave. The leading orderconcentration is uniform across the pipe, so this term is simply:??wr,ave =???? [??(1? ??)??wr]??=??0.Spreading of the front of the slug due to settling occurs when ??wr,avedecreases with ??0, since the lower concentrations are then advected fur-thest. This situation is usually the case, since the settling velocity de-creases with ??0. Similarly, at the rear of a slug we may expect to seesome sharpening of the front.We can estimate the difference in speeds of settling at the front of theslug, by??wr,ave(??0 = 0)? ??wr,ave(??in) ? ??wr,ave(??0 = 0) = ??wr(??0 = 0),402.2. Particle dispersioni.e. this assumes that the settling speed at the input concentration ??inis relatively small. Consequently, a length estimate, for the effect ofspreading due to settling isLmix1 =L?mix1L?=1L?[L?U?0??wr(??0 = 0)]=??wr(??0 = 0)U?0, (2.92)which we have expressed as a ratio of the length of the pipe.(ii) We can also estimate the net effect of diffusivity on spreading the slug.From the Taylor dispersivity DT and the range of ??0 pumped into thewell, we can estimate a representative dispersivity DT,0, (e.g. either themaximum or an average value). A diffusive length-scale is given byL?mix2 =?8D?d,aveL?U?0[1 + Pe2DT,0],or a ratio of the pipe length:Lmix2 =L?mix2L?=?8D?d,aveL?U?0[1 + Pe2DT,0] = 2?1/2?1Pe + PeDT,0.(2.93)In Figs. 2.4-2.5 we have explored the variation in ??wr/U?0, which characterizesLmix1. Although generally small the length Lmix1 is independent of the pipelength L?. On the other hand Lmix2 depends on the pipe aspect ratio, (pro-portionate to ?1/2). The two estimates are plotted in Fig. 2.8. Although formany parameter ranges the diffusive effect will dominate, in any sufficientlylong pipe Lmix1 can be significant with respect to Lmix2. This type of estimateis perhaps the easiest to use for design purposes and note that it could be moreappropriate for application to scale with the distance between slugs, rather thethe pipe length, in assessing whether spreading effects are significant.412.2. Particle dispersion5 10 15 20 2500.10.20.30.40.5U0(m/s)?? 0 0.0002 0.0004 0.0006a)Lmix15 10 15 20 2500.10.20.30.40.5U0(m/s)?? 0.035 0.04 0.045 0.05b)Lmix20.2 0.4 0.6 0.8 100.10.20.30.40.5n?? 0 0.0001 0.0002c)Lmix10.2 0.4 0.6 0.8 100.10.20.30.40.5n?? 0.02 0.03 0.04 0.05 0.06 0.07 0.08d)Lmix24 6 8 10 12 14x 10?400.10.20.30.40.5d?p(m)?? 0.0002 0.0006 0.001e)Lmix14 6 8 10 12 14x 10?400.10.20.30.40.5d?p(m)?? 0.033 0.035 0.037f)Lmix2Figure 2.8: Estimates of the lengths that a proppant slug will spread, basedon particle settling (slip), Lmix1, and particle diffusion, Lmix2. Both lengthscales are made dimensionless with the length of the pipe L?. Parameters areK?0 = 0.1Pa.sn, ??Y 0 = 1Pa, s = 2.65, D? = 62mm, ? = 0.001: a) and b)?p = 0.01, n = 0.5; c) and d) U?0 = 15 m/s, ?p = 0.01; e) and f) U?0 = 15 m/s,n = 0.5. 42Chapter 3ResultsIn this chapter, we present results based on the numerical solution of (2.70).Here, we drop subscript 0 for simplicity. In general, we seek for an algorithmto solve the following transport equation that was introduced in chapter 1.????t +?F (??)?z =??z [D(??)????z ] (3.1)The challenge in solving this type of equation is that both advective slip anddiffusive terms are relatively small. The main transport mechanism is clearlyadvection at the mean speed and in capturing this numerically we induceunwanted numerical errors (typically diffusive and/or dispersive). Some careis needed in ensuring that the numerical errors do not mask the physical effects.Therefore, in the first part of this section we outline the numerical algorithmand present test results only. Thereafter, sample physical results are presentedof relevance to slug dispersion.3.1 Algorithm developmentWe deal first separately with the advective component in ?3.1.1 and the dif-fusive component in ?3.1.2, before combining the methods to solve (2.70) in?3.1.3.433.1. Algorithm development3.1.1 Advection equationWe want to solve the following advection problem??t + F (??)z = 0 (3.2)Adopting a finite volume approach implies that we should study numericalmethods of the form??n+1i = ??ni ??t?z (Fni+ 12? F ni? 12 ), (3.3)where ??ni and ??n+1i are the averaged value of the i-th cell at time tn andtn+1 respectively. Here, Fi? 12 is some approximation to the average flux alongz = zi? 12 :F ni+ 12 ?1?t? tn+1tnF (??(zi+ 12 , t)) dt. (3.4)We use Godunov method (see [23]) to obtain an explicit formula for thenumerical flux (3.4) as followF ni+ 12 = F (??ni , ??n+1i ) =???????min??ni <?<??ni+1F (?) ??ni ? ??ni+1,max??ni+1<?<??niF (?) ??ni+1 ? ??ni .(3.5)This formula is valid also for non-convex flux functions. Implementing theGodunov scheme (3.5) is straightforward. The temporal update (3.3) is simpleto compute as it is explicit, but computing the fluxes can be complicated,since an optimization problem has to be solved locally (see above). If the fluxfunction F has a single minimum at the point ? and no local maxima, thenthe formula (3.5) can be simplified toF ni+ 12 = F (??ni , ??n+1i ) = max(F (max(??ni , ?)), F (min(??ni+1, ?))). (3.6)443.1. Algorithm developmentIf the cell averages are selected as a piecewise constant function in (3.3),this is known to result in an overall first order accuracy of the Godunov-type schemes. However, given the cell averages ??ni at time tn, we can insteadreconstruct linear functions in each cell, which leads to a second-order accuratesolution. We denote the piecewise linear function in the cell ith cell as pi(z),which takes the formpi(z) = ??ni + ?ni (z ? zi). (3.7)where ?ni is a parameter that determines the slope in cell ?i. The properchoice of ?ni satisfies the Total Variation Diminishing (TVD) property. Forexample the following so-called minmod limiter is TVD?ni = minmod( ??ni+1 ? ??ni?z ,??ni ? ??ni?1?z), (3.8)where the minmod function is defined asminmod(a, b) =?????sign(a)min(|a|, |b|) sign(a) = sign(b),0 otherwise.(3.9)A wide variety of other slope limiters such as Superbee, Monotonized central-difference (MC), Van Leer and etc, have also been proposed in the literature,(see e.g. [23] for an overview).Despite the second-order accuracy of the piecewise linear reconstruction(3.7), the first-order temporal accuracy leads to an overall first-order accu-racy, which negates the purpose of using a piecewise linear reconstruction.Therefore, instead, we implement a second-order time stepping method. Weuse Runge-Kutta methods that preserve the TVD property. Such methods aretermed Strong Stability Preserving (SSP) Runge-Kutta methods. In particu-453.1. Algorithm developmentlar, we use the following second-order SSP-Runge-Kutta method??? = ??n + ?tL(??n)???? = ??? + ?tL(???)??n+1 = 12(??n + ????), (3.10)where, the operator L denotes any operator acting pointwise on ??, which heretakes the form:L(??(t))i = ?1?z (Fi+ 12 (t)? Fi? 12 (t)). (3.11)In summary, with all the above ingredients in place, we can state thealgorithm for computing with a second-order scheme. Given cell averages ??iat time level tn , we need to perform the following steps:Step 1 (Reconstruction): ??i, reconstruct the averages to obtain the piece-wise linear function (3.7). Any nonoscillatory slope limiter like the min-mod (3.8), Superbee or the MC limiter can be used. Note that we onlyrequire the edge values ??i at each cell.Step 2 (Flux evaluation): Given the edge values ??i at each cell, we plugthese values into the numerical flux (3.5).Step 3 (Time stepping): For second-order schemes, we use the second-order SSP Runge-Kutta method (3.10). As this method consists of 2stages, steps 1 and 2 must be applied to each stage (e.g., ??n and ??? ).The time step ?t in (3.10) is determined by a CFL condition of the formmaxi|?F (??ni )??? |?t?z ? CFL. (3.12)The above second-order high resolution scheme is TVD as all the three ingre-dients are constructed to ensure this property.463.1. Algorithm developmentAs well as the second-order Godunov-type schemes, we have also imple-mented a second-order extension of the staggered Lax-Friedrichs scheme, asintroduced by Nessyahu & Tadmor [29], (commonly called the NT scheme).This scheme is one of the simplest possible examples of second-order cen-tral schemes and slightly faster compared to the second-order Godunov-typeschemes. Generally we have found better results using the NT scheme. As atest, we have used the NT scheme with different limiters to solve the linear(constant speed) advection problem for initial data consisting of a combina-tion of a smooth, squared cosine wave and double step function. Note that forthe fracturing problem the main advective component is a constant velocity,with smaller settling velocity contribution. With the scaling used in deriving(2.70), it is clear that to simulate the physical process we will advect pulsesfor a unit time and distance, hence the relevance of this example. In order toexamine longer time attenuation we impose periodic boundary condition onour unit domain.Figures 3.1 and 3.2 show the numerical solution after 9.5 and 10 time unitsrespectively. We can see that in all cases the results are very robust with aninsignificant amount of numerical diffusion. The exact solution simply advectsthe intial condition (shown) at unit speed, so that Fig. 3.2 in fact comparesthe exact and numerical solutions at t = 10. Of the 4 limiters tested, theSuperbee limiter appears to resolve the discontinuities better.3.1.2 Diffusion equationIn this section, we explain the numerical method for the diffusive component,i.e. we solve the following quasi-linear parabolic equation.??t = [D(??)??z]z. (3.13)There are many methods used for such equations, which are not difficult tosolve. Here, we use the CPC method of [25] which is a simplified version473.1. Algorithm development0 0.5 100.51MinMod0 0.5 100.51vanLeer0 0.5 100.51MC0 0.5 100.51SuperbeeFigure 3.1: Solution of the linear (constant speed) advection problem in aperiodic domain with unit length by the NT scheme at t = 9.5 (i.e. after 9.5time periods), CFL = 0.4,?z = 0.001. Solid line: initial condition; dashedline: numerical solution.0 0.5 100.51MinMod0 0.5 100.51vanLeer0 0.5 100.51MC0 0.5 100.51SuperbeeFigure 3.2: Solution of the linear (constant speed) advection problem in aperiodic domain with unit length by the NT scheme at t = 10 (i.e. after 10time periods), CFL = 0.4,?z = 0.001. Solid line: initial condition; dashedline: numerical solution.483.1. Algorithm developmentof a Newton Predictor-Corrector (NPC) scheme. This method has a second-order convergence in space. Convergence in time direction is slightly less thansecond-order (i.e. O(?t2??), for some small ? > 0), but the method is simpleto implement and stable.The discretization of (3.13) consists of two stages both of which lead to asystem of linear algebraic equation, that can be solved by a simple tri-diagonalmatrix algorithm in favor of a more costly matrix inversion.2?t(???i ? ??ni ) =1(?z)2?z{D(zi+ 12 ,12(??ni + ??ni+1))?z???i}(3.14)1?t(??n+1i ? ??ni ) =12(?z)2?z{D(zi+ 12 ,12(???i + ???i+1))?z(??ni + ??n+1i )}.(3.15)Here ?z and ?z are backward and forward differencing operators respectively.To test this scheme we consider the following three examples.Example 1:??t = [D(??)??z]z, 0 ? z ? 1 0 ? t ? 2D(z, ??) = exp(z)(2z + exp(?z))[1 + (exp(?z) ? 2)/??]. (3.16)Here the exact solution is ?? = (1+exp(t))(2+exp(?z)). We take boundary andinitial conditions consistent with the given exact solution and compare the twosolutions at successive times. Figure 3.3a shows that the CPC method accu-rately approximates a well behaved solution to a parabolic equation with nonconstant diffusivity coefficient. Although the fracturing problem is quasilinear,the nonlinearity merely results in a bounded variation of the axial diffusivityalong the pipe.493.1. Algorithm development0 0.2 0.4 0.6 0.8 100.050.10.150.20.250.30.350.4z??a)t0 0.2 0.4 0.6 0.8 102468101214z??b)tFigure 3.3: Exact (solid line) and numerical solution (dashed line) of diffusionequation for (a) example 1 and (b) example 2 at t = 0, 0.5, 1., 1.5, 2, ?t =0.005 & ?z = 0.005.Example 2:??t = [D(??)??z]z, 0 ? z ? 1, 0 ? t ? 2 D(z, ??) = ??. (3.17)Here the exact solution is ?? = (z + 0.5)2/(18? 3t). This example has a finitetime blow-up and is used to investigate how well the CPC scheme behavesfor a quasilinear equation with a discontinuity in time. Figure 3.3b showsthe results for this case, again showing a good comparison with the analyticalsolution.Example 3:??t = [D(??)??z]z, ? 4 ? z ? 4, 0 ? t ? 20 D(z, ??) = 0.01 (3.18)The initial condition is a unit square wave of length 2 centered at the origin.We take no flux boundary conditions at the end points of our finite domain.On an infinite domain, the exact solution is ?? = 0.5[erf(z+1)/?4Dt? erf(z?503.1. Algorithm development?4 ?2 0 2 400.20.40.60.81z??a)0 5 10 15 2000.511.5tl sb)Figure 3.4: (a) Exact (solid line) and numerical solution (dashed line) of dif-fusion equation for example 3 at different time. (b) (dashed line) shows thecomputed spreading of the front and (solid line) shows the approximationls =?8Dt.1)/?4Dt]. This example mimics the diffusive behavior of the fracturing prob-lem, in a moving frame of reference. The aim however is to understand whatis the proper characteristic length-scale for diffusion.Note that we have approximated the diffusive spreading length-scale usingthe estimate?8Dt, earlier. The coefficient value 8 is derived from the an-alytical solution by requiring that the value of ?? in the domain be & 2% ofthe maximum initial ??. Figure 3.4a shows the accuracy of the CPC schemein approximating the spreading of a pulse. Figure 3.4b compares the estimate?8Dt against a numerical approximation which computes the 2% limit fromthe numerical solution. We can see that the error is reasonable for times ofO(1), which is what is required for the fracturing process. Note that there area number of different error sources here: (i) the coefficient 8 is not precise forthe analytical solution; (ii) effects of the finite domain; (iii) numerical errorsattributable to the scheme used.513.2. Physical results3.1.3 Advection-Diffusion equationLet?s first consider the numerical solution to an advection-diffusion problemof the form??t + F (??)z = D??zz. (3.19)One standard approach for solving (3.19) is to use a fractional-step or operator-splitting method, in which we solve the following two simpler problem??t + F (??)z = 0 (3.20)??t = D??zz (3.21)Algorithmically, we start with initial data ??ni to obtain the intermediate value???i from (3.20), then we take ???i to obtain ??n+1i from (3.21). This scheme isalmost second-order accurate and we refer readers to [23] for a discussion onthe splitting error and order of accuracy of such methods.The splitting step (3.20) allows us to use second-order the method ex-plained in ?3.1.1 for the advection part without change. We also implementthe scheme in ?3.1.2 to take account of the diffusive part. As a test example,in Fig. 3.5 we show results of a constant speed linear advection and constantdiffusion problem, comparing the numerical solution with the exact solution(from an infinite domain). The exact and numerical solutions are in a verygood agreement.3.2 Physical resultsFinally we show some numerical results, solving (2.70) and showing the evo-lution of proppant pulses (slugs) for fracturing parameters. Figure 3.6a showshow an initial slug of proppant advects along the pipe in the absence of dif-fusion. Over the length of pipe considered, the spreading length due to slip523.2. Physical results0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??inFigure 3.5: Exact (solid line) and numerical solution (dashed line) of the linear(constant coefficient) advection-diffusion equation at different times; parame-ters: D = 0.3831, CFL = 0.4, dz = 0.001.of particles is clearly insignificant. We include the diffusion of solid particlesin Fig. 3.6b, which appears to be the dominant spreading process for the pa-rameters considered (a relatively short pipe and low Re). For the same basicparameters, Fig. 3.7b shows the interaction of two initially separate pulses dueto diffusion.533.2. Physical results0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??ina)0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??inb)Figure 3.6: The evolution of a square wave at different time under transportequation (2.70) (a) excluding the diffusion part (b) including the diffusionpart. Here, we have Re = 3100, B = 0.62, F r = 5, s = 2.65, ? = 0.01, ?p =0.01, ??in = 40%, n = 1.0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??inFigure 3.7: The interaction of two square waves transporting with equation(2.70) at different time. Here, we have Re = 3100, B = 0.62, F r = 5, s =2.65, ? = 0.01, ?p = 0.01, ??in = 40%, n = 1.543.2. Physical resultsStill solving (2.70), but moving now to a set of process parameters in arange given in Table 1.1 and Table 1.2, we show the effects of diffusion ona sequence of pulsed proppant slugs in Fig. 3.8. The pattern of pulse length(and interval between pulses) is more characteristic of fracturing process thanthe previous examples. We see that for these particular parameters the pulsesremain quite square as they travel down the pipe. The last image howeverzooms in to a frontal region (here a backwards step) which reveals that thetransitions between pulses are smooth.As second example is shown in Fig. 3.9 for a typical set of fracturing pa-rameters (i.e. the parameters about which we have varied the dimensionlessgroups previously). As pulse duration we have taken ?t = 0.15 and have takenan interval ?t = 0.15 between pulses. Although the pulses remain distinct,there is considerable rounding of the step discontinuities by the time the pulsesattain the bottom of the pipe. Although there is no interaction at present,reducing the time interval between pulses would soon result in interaction.To illustrate that pulse interaction can occur relatively easily, we haverepeated the computation of Fig. 3.9, but with double the frequency of pulses(hence halving both the pulse duration and the time interval between pulses).Figure 3.10 shows the results of this change, which results in a modest amountof interaction by the end of the pipe.The effects of pipe length are shown in Fig. 3.11 which follows the sameparameters as Fig. 3.9, but with a pipe length 3 times as long (? = 0.00033).We have kept the dimensionless frequency of pulses the same (although ar-guably this should also have increased). We see that the pulses are interactingby the end of the pipe and the last image in Fig. 3.11 zooms into the pulsesat the last time of the snapshot. We can see that as well as interacting thepulses are not completely symmetric. This may be due to the settling effectas discussed earlier.553.2. Physical results0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0.926 0.928 0.93 0.932 0.9340.020.20.40.60.80.98z??/??inFigure 3.8: Results of solving (2.70) for a typical set of process parameters anda dimensionless time interval of 0.1 between pulses. The last figure shows thespreading of the last pulse (backward step), amplified in the axial coordinate.Here, we have Re = 4318, B = 0, F r = 17.52, s = 2.65, ? = 2.067? 10?5, ?p =0.009, ??in = 30%, n = 0.5.563.2. Physical results0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0.6 0.7 0.8 0.900.20.40.60.81z??/??inFigure 3.9: Result of a typical example with Re = 31000, B = 0.31, F r =10, s = 2.65, ? = 0.001, ?p = 0.01, ??in = 30%, n = 1 and pulsation timescale0.15. The snapshots are shown with dimensionless time interval of 0.1. Thelast figure zooms at the last pulses and shows there is no interaction betweenpulses.573.2. Physical results0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0.6 0.7 0.8 0.900.20.40.60.81z??/??inFigure 3.10: We double the number of pulses in previous example (see Fig. 3.9)which results in pulse interaction. The snapshots are shown with dimensionlesstime interval of 0.1. The last figure zooms at the last pulses and shows theinteraction between pulses.583.2. Physical results0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0 0.2 0.4 0.6 0.8 100.20.40.60.81z??/??in0.6 0.7 0.8 0.900.20.40.60.81z??/??inFigure 3.11: In this example the length of the pipe is 3 times longer than inexample of Fig. 3.9. The snapshots are shown with dimensionless time intervalof 0.1. The last figure zooms at the last pulses and shows the interactionbetween pulses.59Chapter 4SummaryThe work presented in this thesis has given a simple model to describe disper-sion of a particulate slug in a vertical pipe, when the flow is turbulent and thefluid is shear-thinning with yield stress. This type of flow finds application inrecent developments of hydraulic fracturing techniques.In chapter 1, we give a brief introduction to fracturing flows and associatedliterature, then we continued with a crude consideration of the types of fractur-ing flow that occur within the pipe. This focuses on bulk flow characteristicsand on what happen at a typical particle scale. In general, over the range offully turbulent suspension flows likely to be observed, we can have flows thatare both inertial and non-inertial on the particle scale, depending largely onthe fracking fluid being pumped and on the flow rates. For those flows thatare non-inertial on the micro-scale, an appropriate simplification would be totreat the slurry as a homogeneous single phase mixture. For higher shear flowsand/or with less viscous fracking fluids the phases are only weakly coupled andare inertial at a local scale, as well as being fully turbulent in the bulk sense.For these flows, the suspension is still largely a homogeneous mixture, but amore sophisticated model is appropriate for describing particle-scale inertialeffects. For both cases, we need to also deal with the non-Newtonian characterof the fracking fluid in developing the model.In chapter 2, we considered a model developed from the two-phase govern-ing equations, assuming that solid and fluid phases can be described as twophases of incompressible continua. We simplified the governing equations toobtain a logarithmic approximation to the velocity profile for the mixture ofsolid and fluid. Also, closures for the slip velocity and friction factor were60Chapter 4. Summaryobtained. We applied the method of multiple timescale to derive a transportequation for solid phase. This transport equation is a 1D advection-diffusionmodel for the leading order volume fraction of solid particles. The advectionpart includes the transport of solid phase with mean flow and particle settling.The diffusive part includes both the average diffusivity and Taylor dispersion.We evaluated the dispersion term for for 3 sensible particle diffusivity models.We showed that the Taylor dispersion effect generally dominates diffusivityand that the spreading of a solids front due to dispersion generally dominatesadvective spreading (rarefaction). Based on this model, we developed 2 esti-mates for the ?mixing length? of a slug of proppant slurry pumped in a clearfracking fluid: due to settling of the particles (advection) and diffusion of theparticles. We gave estimates of these mixing lengths for a wide range of flowparameters.In chapter 3, we provided an accurate numerical algorithm for solution ofthe nonlinear 1D advection-diffusion model. Using this we showed how pulsesof proppant may or may not interact for typical process parameters.In general, we can summarize the significant novel contributions of thisthesis as follows.? The method of multiple timescale has been applied to derive a 1D ad-vection/diffusion model for the leading order volume fraction of solidparticles. The model includes both settling and dispersion effects.? A robust numerical scheme has been developed and programmed to solvethe nonlinear advection-diffusion equation for solid transport.? 2 analytical mixing lengths estimates have been introduced (due to ad-vection (settling) and diffusion of particles), in order to estimate thespreading of pulses more directly, i.e. suitable for field application.Recommendations and future directionsThe following areas appeared to show promise for future development.61Chapter 4. Summary? It is of practical interest to study wellbore inclination effects. In dealingwith inclined and horizontal sections of well prior to the fracture, evensmall settling velocities become important. There are numerous criteriafor predicting flow regimes in horizontal slurry transport, e.g. [44]. Al-though the vertical part of the pipe transport may be relatively uninflu-enced by settling for many flows, the horizontal sections need designingto lie above the critical velocity for bed formation.? In using particle diffusivity model of Eskin [17] and Walton [48], weimplemented our own closure for friction factor, which might not beaccurate. More clearly, the estimates of Taylor dispersion are whollyrelate to the mean axial velocity profile, which are in turn coupled to thefriction factor. Essentially, these should be applied together.? Experimental, computational and theoretical techniques should be con-sidered to improve the available closures for the drag coefficient, settlingvelocity, particle diffusivity and the constitutive laws for yield stress fluidsuspensions.62Bibliography[1] Bu?rger, R. Phenomenological foundation and mathematical theory ofsedimentation-consolidation processes Chemical Engineering Journal. 80,177-188, (2000)[2] Mahaut., F. Chateau, X., Coussot, P. and Ovarlez, G. Yieldstress and elastic modulus of suspensions of noncolloidal particles in yieldstress fluids. J. Rheology. 52, 287-313, (2008)[3] Chateau, X., Ovarlez, G., Luu Trung, K Homogenization approachto the behavior of suspensions of noncolloidal particles in yield stress fluids.J. Rheol. 52, 489-506, (2008)[4] P. Coussot and C. Ancey Rheophysical classification of concentratedsuspensions and granular pastes. Phys. Rev. E 59(4), 4445-4457, (1999)[5] Coussot, P., Tocquer, L., Lanos, C., Ovarlez, G. Macroscopic vslocal rheology of yield stress fluids. J. Non-Newt. Fluid Mech. 158, 85-90,(2009)[6] Crowe, C.T., Troutt, T.R. and Chung. J.N. Numerical models fortwo-phase turbulent flows. Ann. Rev. Fluid Mech. 28, 11-43, (1996)[7] Crowe, C.T. 2006 Multiphase Flow Handbook CRC Press[8] Crowe., C.T, Schwarzkopf., J.D., Sommerfeld., M. and Tsuji.Y2012 Multiphase Flows with droplets and particles CRC Press63Bibliography[9] Crowe., C.T On models for turbulence modulation in fluid-particle flowsInternational Journal of Multiphase Flow 26, 719-727, (2000)[10] Drew, D. A. Mathematical modeling of two-phase flow. Ann. Rev. FluidMech. 15, 261-291, (1983)[11] dH?uteau, E., Gillard, M., Miller, M., Pena, A., Johnson, J.,Turner, M., Medvedev, O., Rhein, T., & Willberg, D. Open-Channel Fracturing A Fast Track to Production Oilfield Review Autumn,pp. 4-17, (2011)[12] Dougherty, T.J. 1959 Some problems in the theory of colloids. PhD Thesis,Case Institute of Technology.[13] Drew, D. A. Mathematical modeling of two-phase flow. Ann. Rev. FluidMech. 15, 261-291, (1983)[14] Eskin, D., Leonenko, Y., Vinogradov, O. On a turbulence model for slurryflow in pipelines. Chem. Eng. Sci. 59, 557-565, (2004)[15] Eskin, D., An engineering model of solids diffusivity in hydraulic convey-ing Powder Technology 159, 78-86, (2005)[16] Eskin, D., Miller, M.J., A model of non-Newtonian slurry flow in a frac-ture Powder Technology 182, 313-322, (2008)[17] Eskin, D., A simple model of particle diffusivity in horizontal hydrotrans-port pipelines Chem. Eng. Sci. 82, 84-94, (2012)[18] Frigaard, I.A. Dispersion of solids in fracturing flows: strategic model-ing review. Confidential consulting report prepared for Schlumberger NTC,April 2013.[19] Hammond, P.S. Settling and slumping in a Newtonian slurry, and im-plications for proppant placement during hydraulic fracturing of gas wells.Chem. Eng. Sci. 50, pp. 3247-3260, (1995)64Bibliography[20] Ishii, M. and Hibiki.T 2011 Thermo-fluid dynamics of two-phase flowSpringer.[21] Krieger, I.M., Rheology of monodisperse lattices. Adv. Coll. Int. Sci. 3,111-136, (1972)[22] Lakhtychkin, A., Vinogradov, O., Eskin, D. Modeling of Place-ment of Immiscible Fluids of Different Rheology into a Hydraulic FractureInd. Eng. Chem. Res. 50, pp. 5774-5782, (2011)[23] Leveque, R.J 2002 Finite volume methods for hyperbolic problems Cam-bridge texts in applied mathematics.[24] Chateau, X., Ovarlez, G. and Trung, K.L. Homogenization ap-proach to the behavior of suspensions of noncolloidal particles in yield stressfluids. J. Rheology. 52, 489-506, (2008)[25] Meek, P.C., Norbury, J. Two-stage, two-model finite differenceschemes for non-linear parabolic equations. IMA. J.Num.Analysis 2, 335-356, (1982)[26] Mahaut, F., Chateau, X., Coussot, P., Ovarlez, G. Yield stressand elastic modulus of suspensions of noncolloidal particles in yield stressfluids. J. Rheol. 52, 287-313, (2008)[27] Maron, SH, Pierce, PE, 1956 J Colloid Sci 11, 80-95[28] Mobbs A.T., Hammond, P.S. Computer Simulations of ProppantTransport in a Hydraulic Fracture SPE Production & Facilities, May 2001,pp. 112-121, (2001); SPE paper 69212.[29] Nessyahu, H., Tadmor, E. Non-oscillatory central differencing for hy-perbolic conservation laws. J. Comput.Phys. 87, 408, (1990)[30] Oroskar, A.R., Turian, R.M. The critical velocity in Pipeline flowof slurries. AIChE. 26, 550, (1980)65Bibliography[31] Ovarlez, G., Bertrand, F., Rodts, S. Local determination of theconstitutive law of a dense suspension of noncolloidal particles through MRI.J. Rheol. 50, 259-292, (2006)[32] Ovarlez, G., Bertrand, F., Coussot, P., Chateau, X. Shear-induced sedimentation in yield stress fluids. J. Non-Newt. Fluid Mech. 177-178, 19-28, (2012)[33] Pearson, J.R.A. On suspension transport in a fracture: framework fora global model. J. Non-Newt. Fluid Mech. 54, pp. 503-513, (1994)[34] Quemada, D. (1982) Lecture Notes in Physics, 164: Casas-Vazquez J,Lebon J (eds), Stability of Thermodynamic Systems, pp 210-247, Springer,Berlin[35] M.C. Roco and C.A. Schook Modelling of slurry flow: the effect ofparticle size. Can. J. Chem. Eng., 61, pp. 494?503, (1983)[36] J.F. Richardson and W.N. Zaki Sedimentation and fluidisation:Part 1. Trans. Inst. Chem. Eng., 32, pp. 35-53, (1954)[37] M.C. Roco and C.A. Schook A modelling for turbulent slurry flow.J. Pipelines, 4, pp. 3?13, (1984)[38] Rudman., M. Blackburn, H.M., Graham, L.J.W., Pullum, L.Turbulent pipe flow of shear-thinning fluids J. Non-Newt. Fluid Mech.118(1), 33-48, (2004)[39] Samir, M. Novel Technique to Increase Production from Tight Reservoirsusing Channel Fracturing Technique for the first time in the Middle East andNorth Africa SPE Sahara Oil & Gas Company, April 2013, pp. 1-10; SPEpaper 164655.66Bibliography[40] Samir, M., Kamal, M., Mathur, A., Semary, M., Yosry, M.,Bernechea, J.M., Kamar, A. & Mokhtar, A. Production Enhance-ment in the Egyptian Western Desert using the Channel Fracturing Tech-nique A Field Case Study SPE Sahara Oil & Gas Company, April 2013,pp. 1-10; SPE paper 164709.[41] Tabuteau H, Coussot P, de Bruyn JR. Drag force on a sphere in steadymotion through a yield-stress fluid. J. Rheol. 51:125-37, (2007)[42] Taylor, G. The dispersion of matter in turbulent flow through a pipe.Proceedings of the Royal Society of London. Series A. Mathematical andPhysical Sciences 223 (1155), 446-468, (1954)[43] Tennekes, H. and Lumley.J.L 1972 A First Course in TurbulenceMIT press.[44] Turian., R.M., Yuan, T.-F. Flow of slurries in pipelines AIChE Jour-nal, 23(3), 232-243, (1977)[45] Truesdell, C. 1984 Rational thermodynamics Springer-Verlag, NewYork-Berlin-Heidelberg.[46] Valco, P. and Economides.M. J 1995 Hydraulic Fracture MechanicsJOHN WILEY & SONS.[47] Vu, T.S., Ovarlez, G., Chateau, X. Macroscopic behavior of bidis-perse suspensions of noncolloidal particles in yield stress fluids. J. Rheol.54, 815?833, (2010)[48] Walton I. C. Eddy diffusivity of solid particles in a turbulent liquidflow in a horizantal pipe. AIChE. 41(7), 1815-1820, (1995)[49] Wierenga, A.M., Philips, A.P., Low-shear viscosity of isotropic disper-sions of (Brownian) rods and fibres; a review of theory and experiments.Coll. Surf. A, 137, 355-372, (1998)67
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Transport and dispersion of particles in visco-plastic...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Transport and dispersion of particles in visco-plastic fluids Hormozi, Sarah 2013
pdf
Page Metadata
Item Metadata
Title | Transport and dispersion of particles in visco-plastic fluids |
Creator |
Hormozi, Sarah |
Publisher | University of British Columbia |
Date Issued | 2013 |
Description | This thesis focuses on development of a model to predict “spreading” of the solids (i.e. proppant) fraction during the fracturing operation. We develop a 1D model that allows us to estimate dispersion of solid particles along a vertical pipe in a fully turbulent flow of a shear thinning yield stress fluid (i.e., visco-plastic fluid), as well as slip relative to the mean flow. In dimensionless form, this results in a quasilinear advection-diffusion equation. Advection by the mean flow, particle settling relative to the mean, in the direction of gravity, turbulent particle dispersivity and Taylor dispersion are the 4 main transport phenomena modelled in the 1D model. We provide a simple analysis of the 1D model, suitable for spreadsheet-type field design purposes, in which we estimate “mixing lengths” due to both settling and dispersion. Secondly, we provide an accurate numerical algorithm for solution of the 1D model and show how pulses of proppant (i.e. slugs) may or may not interact for typical process parameters. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2013-09-03 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0074242 |
URI | http://hdl.handle.net/2429/45002 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 2013-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
Aggregated Source Repository | DSpace |
Download
- Media
- 24-ubc_2013_fall_hormozi_sarah.pdf [ 1.7MB ]
- Metadata
- JSON: 24-1.0074242.json
- JSON-LD: 24-1.0074242-ld.json
- RDF/XML (Pretty): 24-1.0074242-rdf.xml
- RDF/JSON: 24-1.0074242-rdf.json
- Turtle: 24-1.0074242-turtle.txt
- N-Triples: 24-1.0074242-rdf-ntriples.txt
- Original Record: 24-1.0074242-source.json
- Full Text
- 24-1.0074242-fulltext.txt
- Citation
- 24-1.0074242.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0074242/manifest