Journal of Petroleum Science and Engineering 55 (2007) 301 – 316 www.elsevier.com/locate/petrol Averaging model for cuttings transport in horizontal wellbores Gilberto Espinosa-Paredes a , Rubén Salazar-Mendoza b , Octavio Cazarez-Candia b,⁎ a b Departamento de Ingeniería de Procesos e Hidráulica, Universidad Autónoma Metropolitana-Iztapalapa, Av. San Rafael Atlixco 186 Col. Vicentina, México, D.F., C.P. 09340, México Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas Norte No. 152, México D.F., C.P. 07730, México Received 31 May 2005; received in revised form 22 March 2006; accepted 28 March 2006 Abstract This paper presents a theoretical analysis of the problem of cuttings transport for a two-region system composed of a fluid bed (ω-region) and a stationary bed of drill cuttings, which is considered as a porous medium (η-region) in the two-phase system. The ω-region is made up of a solid phase (σ-phase) dispersed in a continuos fluid phase (β-phase), while the η-region consists of a stationary solid phase (σ-phase) and a fluid phase (β-phase). The volume averaging method was applied in this study. Volumeaveraged transport equations were derived for both the fluid bed and the porous medium regions. These equations are based on the non-local form of the volume-averaged momentum transport equation that is valid within the bounded region. Outside this region, the non-local form reduces to the classic volume-averaged transport equation. From these equations, a one-equation model was obtained, and the constraints that the one-equation model must satisfy were applied. For estimating the averaged pressure drop and the averaged velocity, the one-equation model was solved numerically by using the finite-difference technique in the implicit scheme. Numerical results are in agreement with experimental data and theoretical results reported in the literature. © 2006 Published by Elsevier B.V. Keywords: Two-phase; Cuttings transport; Solid–liquid flow; Modeling; Volume averaging 1. Introduction Due to the presence of two phases (solid and liquid) where the solid particles tend to settle at the bottom of the pipe (Doron and Barnea, 1993), the hydraulic transport of solid particles in horizontal pipes is a very complex physical phenomenon. Such phenomenon is relevant in several areas, such as the chemical, mining and oil industries. In the oil industry, horizontal drilling is used to exploit reservoirs exhibiting thin pay zones, to solve the ⁎ Corresponding author. Tel.: +52 55 91758294; fax: +52 55 51198423. E-mail address: ocazarez@imp.mx (O. Cazarez-Candia). 0920-4105/$ - see front matter © 2006 Published by Elsevier B.V. doi:10.1016/j.petrol.2006.03.027 problems related to water and gas conning, to obtain greater drainage area, and to maximize the productive potential in naturally fractured reservoirs. However, a major deterrent in horizontal drilling is the reduction in performance of the transport of solid rocks fragments called cuttings transport (Cho et al., 2000). Therefore, numerous mathematical and empirical models for the prediction of cuttings transport in horizontal and directional wells have been developed by several researchers (Azar and Sanchez, 1997; Cho et al., 2000, 2002; Doron et al., 1997; Kamp and Rivero, 1999; Leising and Walton, 1998; Li and Walker, 1999; Martins et al., 1996; Nguyen and Rahman, 1996, 1998; Pilehvari et al., 1996; Sanchez et al., 1999; Santana et al., 1998). 302 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 A detailed review of published experimental data reveals that the cuttings transport characteristics change with an increase in wellbore angle. Leising and Walton (1998), Peden et al. (1990), and Sifferman and Becker (1992) reported that, under a certain range of well deviation, the cuttings bed in annuli is unstable. Ford et al. (1990) and Tomren et al. (1986) carried out experimental work on cuttings transport in inclined wellbores and observed the existence of different layers that might occur during the mud flow and cuttings in a wellbore: a stationary bed, a sliding bed, and a heterogeneous suspension or clear mud. Later, Li and Walker (1999) presented results of a study of 600 experimental tests. From these, they developed a computer program for predicting cuttings transport in a multiphase gas–liquid-cuttings system. The sensitivity of the cuttings bed height with respect to the liquid/gas volume flow rate ratio, in situ liquid velocity, rate of penetration, inclination angle, and circulation fluid properties were evaluated. Gavignet and Sobey (1989) used the previous experimental observations as a basis for developing a two-layer cuttings transport model. They assumed that the cuttings fall to the lower part of the inclined wellbore and form a bed that slides up the annulus, and that above this bed, a second layer of pure mud exists. They considered the tubing eccentricity in their geometrical calculations of wetted perimeters and calculated an apparent viscosity for non-Newtonian muds using a rheogram formulated in polynomial form. Their work was extended by Sharma (1990) who separated the particle layer into two separate layers: a stationary bed and a sliding bed, or a bed sliding up inside the annulus on top of a bed sliding down at the bottom. His results agreed very well with experimental data, but few details about the model closure relationships were given. On the other hand, Martins and Santana (1992) presented a two-layer model that is more versatile than Gavignet and Sobey's model because it allows the particles to be in suspension in the upper layer. In this layer, they calculated the mean particle concentration from a concentration profile that they obtained by solving a diffusion model. This approach was based on an earlier work on slurry transport carried out by Doron et al. (1987). Martins et al. (1996) extended the model of Martins and Santana (1992) to flows of non-Newtonian fluids by including an apparent viscosity that is calculated from a yield-power law model. They also presented experimental data on erosion of a sandstone bed by polymeric mud in a horizontal wellbore. The model and experimental data were used to assess the validity of the interfacial friction law, the internal angle of friction, and the dry dynamic friction factor. Santana et al. (1998) addressed similar considerations with regard to interfacial friction factors. Some researchers have proposed three-layer models. For example, Nguyen and Rahman (1996) presented a three-layer model that was organized differently from the model presented by Sharma (1990). They distinguished a layer of cuttings of uniform concentration at the bottom of the wellbore that moves with uniform velocity. On top of that layer, a dispersed cuttings layer can be found with a velocity gradient that follows a distribution proposed by Wilson (1987). The third and upper layer consists of clear mud. A sensitivity study of the model on some of its parameters was provided, but the paper did not specify the boundary conditions of each flow mode nor provided quantitative comparison with experimental data. Furthermore, this model did not consider the rheology of the drilling fluid or the sphericity of the cuttings. Doron and Barnea (1993) also worked on a threelayer model, which was obtained from their original work (Doron et al., 1987) on a two-layer model in horizontal pipelines. In their model, they assumed a stationary and a sliding bed of cuttings as well as a heterogeneous suspension. Later, Doron et al. (1997) included gravity acceleration so that the model could be applied to pipelines with a small inclination from the horizontal (less than 10° approximately). Doron et al. (1997) also verified the validity of their model by comparison with experimental data and obtained a good agreement for the dependence of the pressure gradient on total flow rate and on mixture concentration. However, the application of this model has some limitations because it does not consider flow in the annulus, the rheology of the carrier fluid and the rolling/ lifting mechanism for solid transport. In their works, Cho et al. (2000, 2002) extended the mathematical model of Doron and Barnea (1993) and Doron et al. (1997) to include flow in the annulus. They also considered additionally the rheology of the drilling fluid, the drill cuttings shape/concentration and the wellbore geometry with eccentricity of the coiled tubing. On the other hand, Kamp and Rivero (1999) presented a review on the state of the art of mechanistic cuttings transport modeling in highly inclined wellbores. They developed a mechanistic cuttings transport model that is able to predict cuttings bed build-up during drilling; however, their model over-predicts cuttings transport at given mud flow rate (631– 3155 × 10− 6 m3 s− 1). On the basis of the information above, it is fair to say that the common problems associated with most of the G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 existing cuttings transport models relate to inaccurate predictions when compared with experimental data or with in situ drilling results, as well as discrepancies among the models (Azar and Sanchez, 1997). There seem to be two main reasons for these problems. Firstly, the models attempt to cover too wide a range of conditions (from vertical to horizontal). Secondly, the models include too many assumptions or neglect certain observed phenomena (Cho et al., 2002). The aim of this paper is to derive a mathematical model of cuttings transport in horizontal wellbores using the concept of a two-layer solid–liquid flow and the method of volume averaging to predict the flow performance and to evaluate the effects of some important parameters which affect the mechanics of cuttings transport during horizontal well drilling. In order to accomplish the objective of this study, the model is developed using the method of volume averaging (Whitaker, 1999) which is a technique used to rigorously derive transport equations for multiphase systems and one of the main approaches in two-phase flow modeling (Espinosa-Paredes et al., 2002). The model was solved with a numerical approach and two cases were analyzed: (1) fully suspended flow and (2) flow with a stationary bed. The numerical results were compared with experimental data from Doron et al. (1987) and experimental data and theoretical results from Doron and Barnea (1993). 2. Point equations The system under consideration is illustrated in Fig. 1, where the fluid bed system is identified as the ω-region and the porous medium as the η-region. An exploded view of the ω-region that is made up of the solid phase (σ-phase) dispersed in a continuous fluid phase (β-phase) is also shown in Fig. 1. Additionally, an exploded view of the η-region that consists of a stationary solid phase (σphase) and the fluid phase (β-phase) is shown there. Note that the β-phase is flowing in both ω and η regions. The governing point equations, boundary (BC) and initial (IC) conditions that describe the process of momentum transfer in both ω- and η-regions are given by (i) ω-region jd vb ¼ 0 in the b phase  
 Avb qb þ jd vb vb ¼ jpb þ jd T b þ qb g At ð1Þ jd vr ¼ 0 ð3Þ in the b phase in the r phase ð2Þ qr   Avr þ jd ðvr vr Þ ¼ At 303 jpr þ jd Tr þ qr g phase ð4Þ I:C:1 vb ðr; t ¼ 0Þ ¼ f ðrÞ ð5Þ I:C:2 vr ðr; t ¼ 0Þ ¼ gðrÞ ð6Þ in the r B:C:1 vb ¼ vr at the b B:C:2 vb ¼ 0 y¼h B:C:3 ðT b at the b r interface ð7Þ ð8Þ pb IÞd nbr þ ðT r pr IÞd nrb ¼ 0 ð9Þ r interface (ii) η-region jd vb ¼ 0 0¼ in the b phase ð10Þ jpb þ qb g þ lb j2 vb in the b ð11Þ phase B:C:4 vb ¼ 0 at the b B:C:5 vb ¼ 0 y¼ B:C:6 jd hvb i ¼ 0 r interface ð13Þ H y¼ ð12Þ H ð14Þ where v, ρ, T are local variables representing the velocity vector, the density and the total stress tensor (laminar and turbulent), respectively, p is the pressure and g is the gravity acceleration vector. In the jump condition given by Eq. (9), nβσ is the unit vector normal to the interface pointing out of the β-phase. In the solid phase, the term pσ in Eq. (9) represents the solid phase pressure, which is due to collisional and kinetic effects (Huilin and Gidaspow, 2003). The inertial effects (ρβ∇· vβvβ) in the η-region were negligible, while these effects were included in the ωregion. Also, in the η-region, Newtonian flow is postulated. The boundary condition at y = − H has been expressed by Eq. (14) in a form that is suitable for use with Darcy's law; thus, we have used 〈vβ〉 to represent the superficial volume-averaged velocity. The boundary conditions given by Eqs. (8) and (14) are an indication of the mismatch of the length scales that one often encounters in transport problems that involve porous media. The point boundary condition given by Eq. (8) is based on the idea that the point velocity is continuous, while the volume-averaged boundary condition given by Eq. (14) is an approximation based on the idea that the interface at y = − H is impermeable. 304 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 Fig. 1. Cutting transport system and averaging volume. As stated previously, two homogeneous regions are distinguished: a homogeneous ω-region and a homogeneous η-region both constituted by the σ-phase and the β-phase (Fig. 1). We denominate homogeneous regions those portions of the system that are not influenced by the rapid changes that occur in the boundary between regions. Then, to study the process that takes place in two-phase flow, volume-averaged transport equations that are valid within both homogeneous regions must be developed. Thus, the averaging model for the ω-region that describes cuttings transport at the macroscopic level is developed in this paper. In order to describe the process of cuttings transport illustrated in Fig. 1, the volume-averaged form of Eqs. (1)–(4) are developed (ω-region). Development of such mathematical model for this system is relatively straightforward when classic length-scales constraints are satisfied (Carbonell and Whitaker, 1984; Zanotti and Carbonell, 1984); however, difficulties arise in the neighborhood of the ω–η boundary where there exist rapid changes in the liquid volume fraction and the length-scale constraints fail. 3. Averaging equations Using the volume averaging theorems, given in Appendix A, on Eq. (1) the volume averaged of the continuity equation for the β-phase is The procedure leading to the β-phase continuity equation can be repeated for the σ-phase beginning with Eq. (3), and the result is given by Aer þ hvr ir d jer þ er jd hvr ir ¼ 0: At Turning our attention to the momentum equation given by Eq. (2), also using the volume averaging theorems given in Appendix A, the superficial-averaged momentum transport equation is given by  A qb eb hvb ib þ qb jd hvb vb i ¼ jhpb i þ jd hT b i |fflffl{zfflffl} |fflfflfflffl{zfflfflfflffl} At |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} pressure stress tensor accumulation Z convection Z 1 1 nbr d T b dA nbr pb dA þ eb qb g þ ð17Þ _ Abr |fflffl{zfflffl} _ Abr gravity |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} interfacial force where ρβ is considered as a constant within the averaging volume, and εβ = 〈1〉. It is important to note that no length-scale constraints have been imposed on the volume-averaged transport equation. The absence of any length-scale constraint simply means that Eq. (17) is also valid in the boundary between the ω- and η-regions. Eq. (17) can be written in compact form as (Salazar-Mendoza, 2004) qb Aeb þ hvb ib d jeb þ eb jd hvb ib ¼ 0: At  b  A eb hvb ib þ qb jd eb hvb i hvb ib þ qb jd hve b ve b i At   þ qb jd hvb vb iexc ¼ ð15Þ ð16Þ þ eb qb g þ M br eb jhpb ib þ eb jd eb 1 hT b i ð18Þ G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 305 in which the vector Mβσ is defined as Z h i 1 e b jxþy dA M br ¼ nbr d Ip eb jxþy þ T _ Abr Z h   1 nbr d I hpb ib jxþy hpb ib jx þ _ Abr  i þ hT b ib jxþy hT b ib jx dA ð19Þ and represents the interfacial force per unit volume applied on phase β, whereas the excess convective term is defined as hvb vb iexc ¼ hvb vb i eb hvb ib hvb ib hve b ve b i: ð20Þ In Eqs. (18)–(20), ψ̃ (were ψ is some function) represents the spatial deviations around averaged values of the local variables. In Eq. (19), the variables with subscripts x + y are evaluated at points within the averaging volume not located at the centroid x (Fig. 2). The procedure leading to the β-phase momentum equation can be repeated for the σ-phase beginning with Eq. (4), and the result is given by A qr ðer hvr ir Þ þ qr jd ðer hvr ir hvr ir Þ þ qr jd hve r ve r i At
 þ qr jd hvr vr iexc ¼ er jhpr ir þ er jd er 1 hTr i ð21Þ þ er qr g þ M rb where the vector Mσβ represents the interfacial force per unit volume applied on phase σ, and 〈vσvσ〉exc has an expression similar to Eq. (20). Fig. 2. Averaging volume. where Aσ is the projected area of a particle and the vector vr is the relative velocity given by vr ¼ hvr ir hvb ib : ð24Þ In Newton's regime, the drag coefficient is given by (Ishii and Mishima, 1984) CD ¼ 0:45 1 þ 17:67f 6=7 18:67f 2 ð25Þ where f ¼ 1  er Þ1=2 1 er  0:62 1:55 : ð26Þ The virtual mass force is given by FV ¼ CV qb av in the homogeneous x When certain length-scale constraints are satisfied, i.e. lσ, lβ ≪ ro ≪ L, the force Mσβ has an especially simple form in a homogeneous fluid bed (Ishii and Mishima, 1984): where av is the virtual mass acceleration and is given (Drew et al., 1979) by M rb 1 ¼ ðer F D þ er F V Þ Vrx ð22Þ where FD, FV and Vσω are the drag force, the virtual mass force and the volume of the σ-phase in the ωregion. The drag force acting on the particle under steady-state conditions can be written in terms of the drag coefficient CD: 1 CD qb vr jvr jAr 2 in the homogeneous x FD ¼ region ð23Þ av ¼ db hvb ib dt dr hvr ir dt region ð27Þ 3.1. Closure relationships ð28Þ here dα/dt = ∂/∂t + 〈vα〉α·∇ for α = β, σ. The virtual mass coefficient, also known as the volume associated with the induced mass Cv for single particle in an infinite medium can be obtained from potential theory. In order to apply this theory, it is necessary impose the following length-scale restrictions (Espinosa-Paredes, 2001): lδ ≪ lβ,lσ ≪ ro.lδ is the characteristic length of the boundary layer thickness, i.e., the interfacial friction effects are taken as concentrated in a thin boundary layer region around 306 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 the solid cuttings, where viscous effects are important. The flow outside the boundary layer is considered ideal. Then, for a spherical particle, Cv = 1/2, this result was obtained with a concentric cell model (EspinosaParedes, 2001). In order to obtain the vector Mβσ that represents the interfacial force per unit volume applied on phase β, it is considered that accumulation effects are insignificant. Then, the momentum interfacial transfer between the solid and liquid phases is given by Eq. (22) for each phase and can be rewritten as M br ¼ M rb on Abr of the homogeneous x region: ð29Þ The intrinsic average of the dyad of the velocity deviations for the β-phase (〈ṽβṽβ〉β) can be obtained for a homogeneous fluid bed when the length-scale constraints lσ, lβ ≪ ro ≪ L are satisfied and is given by hve b ve b ib ¼ Ker m2r for the homogeneous x region ð30Þ where, for a concentric cell, K represents the stress second-order tensor due to spatial deviation, which is given by (Wallis, 1989) K¼ 1 3 3 ex ex þ ey ey þ ez ez : 5 4 4 ð31Þ The intrinsic average of the dyad of the velocity deviations for the σ-phase (〈ṽσṽσ〉σ) may be null if is considered that in a homogeneous fluid bed all cuttings displace at the same velocity. Then, hv̂r v̂r ir ¼ 0 for the homogeneous x region: ð32Þ With this result, Eq. (21) is further simplified. The averaging momentum transport equations that describe the transport phenomena in the η-region can be obtained from Eqs. (10) and (11) and is given by Whitaker (1986) as 0¼ ð33Þ b 2 eb jhpb i þ eb qb g þ lb j hvb i lb ðjeb Þd ½jðeb 1 hvb iފ lb Φb It has been shown elsewhere (Whitaker, 1986) that Φβ has an especially simple form in a homogeneous porous medium provided that certain length-scale constraints are satisfied. Whitaker (1986) derived an expression for Φβ which is given by Φb ¼ K b 1 d hvb i for the homogeneous g ð34Þ region ð36Þ which is valid when the following three length-scale constraints are imposed: ro ≪1 Le Lp1 ro ≪1 Le Lm2 lb ≪ro : ð37Þ In these equations, Kβ represents the Darcy's law permeability tensor, Lε is the characteristic length associated with ∇εβ, Lp1 is the characteristic length associated with ∇〈pβ〉β, and Lν2 is the characteristic length associated with ∇2〈vβ〉. When the constraints indicated by Eq. (37) are valid, the second Brinkman correction is negligible, compared with the first Brinkman correction: lb ðjeb Þd ½jðeb 1 hvb iފ ¼ 0: Use of Eqs. (36) and (38) in Eq. (34) leads to  Kb  hvb i ¼ d jhpb ib qb g eb 1 lb j2 hvb i lb in the homogeneous g 3.2. Averaging equations for the η-region jd hvb i ¼ 0 where the first viscous term is known as the first Brinkman correction, the viscous term involving the gradient of the porosity is known as the second Brinkman correction, and Φβη is a vector defined by Z h   1 nd I pbg jxþy hpbg ibg jx lbg Φbg ¼ Vbg Abr  i þ lbg jvbg jxþy jhvbg ibg jx dA: ð35Þ region ð38Þ ð39Þ where Kβ = Kzzezzez. The Blake–Kozeny equation (Bird et al., 2002) is used to express permeability as Kzz ¼ dp2 e3b kð1 eb Þ2 ð40Þ where dp is the effective particle diameter, and λ is a constant which is obtained from experimental tests. Whitaker (1996) reported that λ = 180, whereas Bird et al. (2002) reported that λ = 150. 307 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 The ω subscript indicates the homogeneous region where Eq. (45) is strictly valid. When the following length-scale constraints are imposed 4. One-equation model for the homogeneous ω-region In the concept of a one-equation model, the mixture is considered as a whole, rather than as two separate phases. It is evident that the one-equation model is simple with respect to the two-phase model. The oneequation model consists of a reduced number of the total averaged equations and closure relationships required in the complete formulation. In order to obtain such oneequation model for the region occupied by the fluid bed, the principle of local pressure equilibrium is introduced, defined as r b fpgx ¼ er hpr i þ eb hpb i : ð41Þ The previous equation has been written on the basis of the principle of local pressure equilibrium, which implies that the process can be characterized by a single pressure. Also, in the simplest view of two-phase flow, the phases are considered to be homogeneously mixed and traveling at the same velocity trough the system. Thus, the process can be characterized by a single velocity: fvgx ¼ hvr ir ¼ hvb ib ð42Þ When Eq. (42) is applied in Eqs. (18) and (21), the dispersion terms are zero, i.e., hve b ve b i ¼ 0 hve r ve r i ¼ 0 ð43Þ ð44Þ The interfacial force terms in Eqs. (18) and (21) are equal and opposite, and they cancel if these two transport equations can be added to obtain a oneequation model. With these considerations, the oneequation model can be written as A fvgx þ qx jd ðfvgx fvgx Þ¼ jfpgx þ jd fTgx þ qx g At for the homogeneous x region qx where qx ¼ er qr þ eb qb fTgx ¼ hT r i þ hT b i ð46Þ ð47Þ l r ¼ ro ; l b ¼ ro ; ro2 ¼ L2 ð48Þ then the excess terms 〈vβvβ〉exc and 〈vσvσ〉exc of Eqs. (18) and (21), respectively, are negligible compared with the other terms of these equations. 5. Coupling of the averaging models of the ω and η-regions A two-region model may be used to describe cuttings transport in a wellbore. This is supported by experimental observations on cuttings transport in annuli (Tomren et al., 1986), where stationary and moving cuttings beds were observed below a heterogeneous cuttings suspension or clear mud. In addition, good results for slurry transport in pipes have been obtained with a two-layer model (Doron et al., 1987). However, many of these models have been constructed on an intuitive basis rather than on the basis of a rigorous analysis of the governing point equations and boundary conditions. Intuitive analysis leads to hidden assumptions and unsupported simplifications. Before writing the final equations for the two-region model, it is necessary to make some comments about some equations obtained in previous sections. Eqs. (1)–(4) are the governing local equations which describe the process of momentum transfer in the ωregion for the β and σ phases at a point. On the other hand, Eqs. (15), (16), (18) and (21) are the governing averaged-equations which also describe the process of momentum transfer in the ω-region for the β- and σ phases, but they do it in terms of averaged variables, that is, using a scale of the order of ro. However, the variables in Eq. (45) which represent the one-equation model, are valid for scales larger than ro, i.e., it is possible to add terms resulting from a force balance on each layer. Then, according the force balances made by Doron and Barnea (1993) and Doron et al. (1987), Eqs. (45) and (39) may be rewritten in the following form: qx A fvgx þ qx jd ðfvgx fvgx Þ¼ jfpgx þ jd fTgx þ qx g At fTgxg d nxg fTgxw d nxw þ ð49Þ þ DHxg DHw 0¼ eb jhpb ibg þ eb qb g þ lb j2 hvb ig fTgxg d ngx DHxg lb K bg1 d hvb ig ð50Þ 308 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 where {T}βw is the wall stress in the ω-region, {T}βη is the stress in the inter-region between the ω and ηregions, nωw(=− nwω) is the unit vector normal to the wall pointing out of the ω-region and nηω(=− nωη) is the unit vector normal to the inter-region pointing out of the η-region, as illustrated in Fig. 3; DHw is the hydraulic diameter in the ω-region and DHωη is the hydraulic diameter in the inter-region. The final formulation for the two-region model is therefore represented by a set of four equations: Eqs. (49) and (50), given above, and the equations given by jd fvgx ¼ 0 in the x region ð51Þ jd hvb ig ¼ 0 in the g region ð52Þ This set is complemented by the inter-regional boundary condition, which is given by B:C:1 fpgx ¼ hpb ibg ; at Axg ð53Þ On the other hand, it is noted that in the ω–η boundary, εβ and εσ undergo significant changes over a distance equal to the radius of the averaging volume r0, as illustrated in the Fig. 4. In this figure, δ represents the thickness of the interfacial region where there exist rapid changes in εβ and εσ. Also in the ω–η boundary, Fig. 5 shows the continuity of the velocity of each region. Thus, further contributions will deal with the Fig. 4. Volume fraction variation in the neighborhood of the nonhomogeneous zone. development of the momentum jump conditions between the ω- and η-regions for cutting transport. 6. Results and discussion The mathematical model given by Eqs. (49)–(53), was used to study two cases: (Case I) fully suspended flow and (Case II) flow with a stationary bed. In the first case, all the cross section of the pipe is occupied by a solid–liquid dispersed flow with a low concentration of drilling cuttings (like the ω-region in Fig. 1), then Eq. (49), without the stress in the inter-region, {T}ωη, and Eq. (51) are used to simulate the flow. For the second case, a stationary bed of drilling cuttings (porous medium) is at the bottom of the pipe above which flows a flow similar to the Case I (Fig. 1), then to simulate the total flow Eqs. (49)–(53) are used. The models for Case I and Case II in one-dimensional form are given in Appendix B. In order to solve the mathematical models, a backward finite-difference implicit scheme with a point-distributed grid was applied. The set of equations for mass and momentum conservation in a discretized form can be written in matrix form as Ad x ¼ b Fig. 3. Unit vectors. ð54Þ in which A is the matrix's coefficients, x is the vector of the dependent variables and b is the vector of the known G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 which is due to the prediction was done with a homogeneous model, where the mixture properties are used. The predictions are compared with some experimental data from Doron et al. (1987), which were measured for superficial velocities approximately from 0.2 to 1.8 m/s. For these velocities the flow pattern observed was principally a flow with a moving bed, whereas only the data for superficial velocities higher than 1.6 m/s correspond to the fully suspended flow pattern. Therefore, the predictions are better when superficial velocities tend to high values. However, the predictions tend to move away with respect to experimental data at low superficial velocities. Then, obviously the model given by Eqs. (B-1) and (B-2) is not adequate to simulate flow with a moving bed or a flow with a stationary bed, but it is appropriate for a fully suspended flow prediction. Fig. 5. Continuity of the global spatial average velocity. parameters. The coefficients of the Eq. (54) are defined in Appendix C for both study cases. 6.1. Case 1. Fully suspended flow The simulated physical system is a pipe of 4.13 m long and 0.05 m in diameter. The simulation was performed for a fully suspended flow. This type of flow pattern is presented at high superficial velocities (higher than 1.6 m/s). The liquid phase used in the simulation is water and the solid phase consists of particles with 0.003 m of diameter and a density of 1240 kg/m3 (data given in Doron et al., 1987). The pressure at the entrance of the pipe is 151,988 Pa and the total volume fraction of solid particles varies from 0.042 to 0.155. The aim of the numerical simulations is to obtain the pressure and velocity profiles as a function of time and space. With the idea of comparing the calculated pressure drop with experimental data, the dimensionless gradient pressure (ΔP)u204E/u (expressed in meter of water by meters of pipe, m/m) is calculated with the next equation: ðDPÞ⁎ ¼ ðPin Pout Þ=L qwater gr 309 6.2. Case 2. Flow with a stationary bed Suppose that a fully suspended flow in a horizontal pipe is Case I. If the slurry flow rate is reduced, the solid particles tend to form a moving bed at the bottom of the pipe (due to density that is higher than that of the carrier fluid). Decreasing the flow rate further, the solid particles tend to form a stationary bed instead of a moving bed. In this work, the last flow pattern mentioned above was simulated using the same data that for Case I and a value for the total volumetric fraction of solid particles of 0.097. With the model given by Eqs. (B-9) (B-10) (B-11) (B-12), and values of 0.52 and 0.65 for εση, Fig. 7 was obtained. This figure shows the relation H/D as a ð55Þ where Pin and Pout are the pressure at the entrance of the pipe and the calculated pressure at the exit of the pipe, which length is L; ρwater is the water density, and gr is the gravity in the radial direction. In Fig. 6, the dimensionless pressure gradient vs. velocity is presented. It can be observed that the prediction for various total volume fractions of drilling cuttings is similar that the results for a flow of liquid, Fig. 6. Dimensionless pressure gradient vs. velocity for a solid–liquid dispersed flow (ρσ = 1240 kg/m3, ρβ = 1000 kg/m3, dp = 3 mm, L = 4.135 m and D = 50 mm). 310 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 function of εσ and εση, obtaining the flow pattern transition for a solid–liquid two-phase flow. From Fig. 7, it can be observed that (1) the point (0,0) represents one phase flow, (2) the points (εσ,0) represent fully suspended flow, (3) the points (0.52,1) and (0.65,1) represent flow through a porous medium (stationary bed) and (4) with the other points; it is possible to predict if the flow is with a stationary bed or with a moving bed and if the upper region is a fully suspended flow or only a flow of liquid. In this case, the dimensionless gradient was calculated by ⁎ þ eg ðDPÞ⁎ ðDPÞ⁎ ¼ ex ðDPÞx g ð56Þ where (ΔP)ω⁎, (ΔP)η⁎, εω and εη are the dimensionless pressure gradients and the volume fractions for the ω and η regions, respectively. In Fig. 8, the transition between flow with a moving bed (right side of the continuous line) and flow with a stationary bed (left side of the continuous line) are shown. In this figure, we can observe that numerical results for (ΔP)⁎ are in agreement with the profile of experimental data and the Doron and Barnea model in the velocity range from 0.6 to 0.95 m/s. The maximum error calculated between the numerical and experimental data was 5%. With the idea of evaluating the behavior of εσω as a function of: εσ, εση and H/D, the next equation was used
 ð57Þ erx ¼ er erg eg =ex : In Fig. 9, all the points on the εσω-axis (i.e., H/D = 0) represent fully suspended flow and all the points on the Fig. 8. Dimensionless pressure gradient vs. velocity for a flow with a stationary bed. Comparison of two-region model against experimental data and three-layer model from Doron and Barnea (1993) (ρ σ = 1240 kg/m3 , ρ β = 1000 kg/m3 , d p =3 mm, L = 4.135 m, D = 50 mm and εσ = 0.097). H/D-axis (i.e., εσω = 0) represent flow through a porous medium (stationary bed) with a flow of only liquid phase at the top of the pipe. The other points on solid and dashed lines represent flow with a stationary bed (like in Fig. 1). Finally, in Fig. 10, the dimensionless pressure gradient as a function of volume averaged velocity and the relation H/D is presented. It can be observed that, for a constant velocity, if the relation H/D grows up, then also, the dimensionless pressure gradient grows up. In practice and experimentally, it is very difficult to maintain the relation H/D constant; however, with Fig. 10, it is possible to predict how the behavior of the dimensionless pressure gradient is as a function of volume-averaged velocity. 7. Conclusions Fig. 7. Flow pattern transition for a two-region model as a function of the total volume fraction of drilling cuttings and the maximum packing at the bottom of the pipe. In this work, the process of cutting transport for a system composed by two regions has been described: the ω-region, which is a fluid bed, and the η-region, which is a porous medium system. A rigorous mathematical model in which each variable is precisely defined has been derived. In order to accomplish this, volume-averaged transport equations were derived for the fluid bed and the porous medium regions. From these equations, a one-equation model was obtained and the constraints that the model must satisfy are identified. Specifically, the coupling conditions between the homogeneous ω-region (σ-phase and β-phase) and the homogeneous η-region (σ-phase) were identified. G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 311 flow with a stationary bed; (2) the model for a flow with a stationary bed [Eqs. (49)–(53)] allows to obtain the flow pattern transition, as a function of total volume fraction, which depends on the maximum packing at the bottom of the pipe, (3) the results for the dimensionless pressure gradient, obtained for a flow with a stationary bed, are in agreement with experimental data and theoretical results reported in the literature, as is illustrated in Fig. 8. Fig. 9. Behavior of volume fraction of drilling cuttings at the top of the pipe as a function of (1) the total volume fraction of drilling cuttings, (2) the maximum parking at the bottom of the pipe, and (3) the relation H/D. The coupled averaged model of ω and η regions given by Eqs. (49)–(53) allows to obtain average pressure drop and average velocity and two cases were analyzed: (1) fully suspended flow and (2) flow with a stationary bed. The numerical results and their comparison with experimental data and theoretical results show that (1) the model for a fully suspended flow is not adequate for simulating a flow with a moving bed nor a Nomenclature A Cross-sectional area A Matrix's coefficients aν virtual mass acceleration =Aσβ area of the β–σ interface contained within Aβσ the averaging volume Aσ area of the σ phase drag coefficient CD virtual mass coefficient of the concentric cell Cν model dp effective particle diameter hydraulic diameter DH normal unit vector in x direction ex normal unit vector in y direction ey ez normal unit vector in z direction drag force FD virtual mass force Fν f friction factor g acceleration vector of gravity H height of the η-region h height of the ω-region Fig. 10. Comparison of the dimensionless pressure gradient as a function of the averaged superficial velocity {vz} and the relation H/D. 312 I j Kβ Kzz lβ lσ lδ L Lε Lp1 Lv2 Mβσ nβσ nηω nωw pβ pσ p̃ (Δp)⁎ pin pout {p}ω ro r Re S t Tβ Tσ T̃ {T}ω vβ vσ vr ṽβ ṽσ {v}ω ∨ Vβ G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 unit tensor unit vector Darcy's law permeability tensor permeability defined by Eq. (40) characteristic length for the β-phase characteristic length for the σ-phase characteristic length for the boundary layer thickness large length scale characteristic length associated with changes in ∇εβ characteristic length associated with changes in ∇〈pβ〉β characteristic length associated with changes in ∇2〈vβ〉 (=− Mσβ) interfacial force per unit volume applied on the β-phase (=− nσβ) normal unit vector directed from the β-phase towards the σ-phase (=− nωη) normal unit vector directed from the ω-region towards the η-region (=− nwω) is the unit vector normal to the wall pointing out of the ω-region pressure in the β-phase pressure in the σ-phase spatial deviation of the pressure dimensionless pressure gradient pressure at the entrance of the pipe pressure at the exit of the pipe pressure defined by Eq. (41) for the oneequation model radius of the averaging volume position vector Reynolds number wetting perimeter time total stress tensor (laminar and turbulent) associated with the β-phase total stress tensor (laminar and turbulent) associated with the σ-phase spatial deviation of total stress tensor total stress tensor defined by Eq. (47) for the one-equation model velocity vector in the β-phase velocity vector in the σ-phase relative velocity vector spatial deviation velocity in the β-phase spatial deviation velocity in the σ-phase velocity defined by Eq. (42) for the oneequation model averaging volume volume of the β-phase contained in the Vσ w x x y 〈〉 {} 〈 〉k 〈 〉exc k 〈 〉m averaging volume volume of the σ-phase contained in the averaging volume interface velocity position vector locating the centroid of the averaging volume vector of the dependent variables position vector relative to the centroid of the averaging volume denotes superficial average denotes variables in the one-equation model denotes intrinsic average for phases k = β, σ denotes excess terms denotes intrinsic average for phases k = β,σ within of the m = ε, η region. Greek symbols β fluid phase δ boundary layer thickness εβ volume fraction of the β-phase volume fraction of the σ-phase εσ a vector given by Eq. (36) Φβ η porous medium λ constant μβ viscosity in the β-phase density in the β-phase ρβ density in the σ-phase ρσ ρω density defined by Eq. (45), for the oneequation model σ solid phase ω fluid bed system ψβ some function associated with the β-phase Subscripts ω identifies a quantity associated with the ω region η identifies a quantity associated with the η region ωη identifies a quantity associated with the ω–η boundary βη identifies the fluid phase in the η-region βω identifies the fluid phase in the ω-region ση identifies the solid phase in the η-region σω identifies the solid phase in the ω-region Acknowledgements The authors acknowledge the financial support provided by the Competencia de Producción de Hidrocarburos of the Instituto Mexicano del Petróleo 313 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 (IMP) and Consejo Nacional de Ciencia y Tecnologia (CONACyT) of México. Appendix A. Volume averaging The superficial volume average of some function ψβ associated with the β-phase is defined as hwb ijx ¼ 1 _ Z Vb ðx;tÞ
 wb x þ yb ; t dVy ðA 1Þ where Vβ (x,t) is the volume of the β-phase contained within the averaging volume ∨ illustrated in Fig. 2. In this figure, it is indicated that x represents the position vector locating the centroid of the averaging volume, while yβ represents the position vector locating points in the β-phase relative to the centroid. In Eq. (A–1) dVy is used to indicate that the integration is carried out with respect to the components of yβ, and the nomenclature used in Eq. (A–1) clearly indicates that volumeaveraged quantities are associated with the centroid. In Eq. (A–1), ∨ = Vσ + Vβ and is independent of space and time; however, Vσ and Vβ depend on x and t. In order to simplify the notation, the precise nomenclature used in Eq. (A–1) is avoided and the superficial average of ψβ is represented as Z 1 hwb i ¼ w dV ðA 2Þ _ Vb b while the intrinsic average is expressed in the form Z 1 w dV ðA 3Þ hwb ib ¼ _b Vb b Both of these averages are used in the theoretical development of this paper. They are related by hwb i ¼ eb hwb i b ðA 4Þ averaging theorem and general transport theorem given by Gray and Lee (1977): Z 1 nbr wb dA ðA 6Þ hjwb i ¼ hjwb i þ _ Abr Awb At Ahwb i ¼ At 1 _ Z wb wd nbr dA ðA 7Þ Abr where w is the interface velocity. Appendix B. One-dimensional models B.1. Case I: Fully suspended flow In this case, Eqs. (49) and (51) in one-dimensional form can be simplified to obtain Afmz gx ¼0 Az qx A d fszz gxw ¼0 fmz gx þ fpgx þ At dz DHx ðB 1Þ ðB 2Þ To obtain Eqs. (B-1) and (B-2), Eq. (51) was used into Eq. (49), the total stress tensor (∇·{T}ω) was neglected and the stress in the inter-region ({T}ωη) does not exist. In Eqs. (B-1) and (B-2), the macroscopic superficial averaged velocity, {νz}, pressure, {pz}, and density, ρω, can be calculated from Eqs. (42), (41) and (46), respectively, whereas the wall stress, {τzz}ωw, is given by (Doron et al., 1987; Doron and Barnea, 1993) 2 1 fszz gxw ¼ fx qx ðfmz gx 2 ðB 3Þ ðB 4Þ where (Doron et al., 1987) fx ¼ aRex f In Eq. (B-4) α = 0.046; ζ = 0.2 for turbulent flow and α = 16; ζ = 1 for laminar flow, and The liquid volume fraction εβ is defined explicitly as eb ðx; t Þ ¼ Rex ¼ Vb ðx; tÞ _ ðA 5Þ It should be clear that the liquid volume fraction εβ is a function of position and depends of the sampling point located by x. In addition to the definitions given by Eqs. (A–2) (A–3) (A–4) (A–5), it is necessary to use the spatial qx fmz gx DHx lx ðB 5Þ where the viscosity is given by (Ishii and Mishima, 1984)  erx  1:625 lx ¼ lbx 1 ðB 6Þ 0:65 where μβω is the viscosity of the β-phase in the ω-region. 314 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 For this case, the volume fraction of the σ-phase in the ω-region is given by exr ¼ Arx A ðB 7Þ In Eq. (B-7) Aσω is the area occupied by the σ-phase and A is the area of the pipe. The hydraulic diameter can be calculated by DHx 4A ¼ S ðB where 〈νβz〉η is the β-phase velocity in the η-region and the friction factor is calculated by (Televantos et al., 1979) 8Þ 1 pffiffiffiffiffiffiffiffiffi ¼ 2fxg 0 1 dp BDHx 2:51 C ffiC 0:86lnB @ 3:7 þ Re pffiffiffiffiffiffiffiffi 2fxg A x ðB 16Þ where dp is the particle effective diameter. In Eqs. (B-10) and (B-12), DHωη is the hydraulic diameter associated with the ω–η boundary and is given by where S is the wetting perimeter of the pipe. DHxg ¼ 4Ax Sxg for the x DHxg ¼ 4Ag Sxg for the g region ðB 17Þ ðB 18Þ B.2. Case II. Flow with a stationary bed In this case, Eqs. (49) and (51) in one-dimensional form can be written respectively as d m z gx ¼ 0 dz qx ðB A d fszz gxw fszz gxg fmz gx þ fpz gx ¼ At dz DHx DHxg 9Þ In Eq. (B-12), Kβηzz is the permeability given by the Blake–Kozeny equation (Eq. (40)), the volume fraction of the β-phase is eb ¼ 1 ðB er ¼ ebx þ ebg ðB 19Þ ðB 20Þ 10Þ and the volume fraction of the σ-phase is er ¼ 1 d hmbz ig ¼ 0 dz region ðB eb ¼ erx þ erg : 11Þ Appendix C. Coefficient vectors d eb hpbz ibg ¼ dz 1 lbKbgzz hmbz ig þ fszz gxg DHxg ðB 12Þ In this case, {τzz}ωw, fω, and Reω are calculated with Eqs. (B-3) (B-4) (B-5), respectively. The volume fraction of the σ-phase and the hydraulic diameter in the ω-region are given by In this section, the coefficient vectors of the finitedifference equations system (Eq. (54)) corresponding to the analysis of each cases are defined. C.1. Case I. Fully suspended flow A¼ " 1 1 Dt # 1 ; qx Dz ðC 1Þ x¼ # fvz gx jtþDt i ; fpz gx jtþDt i ðC 2Þ 14Þ " The stress in the inter-region is given by (Doron et al., 1987; Doron and Barnea, 1993) 2 fvz gx jitþDt 1 erx ¼ ðer A erg Ag Þ=Ax ðB 13Þ and DHx ¼ 4Ax Sx þ Sxg  1 fszz gxg ¼ fxg qx fmz gx 2 ðB hmbz ig 2 ðB 15Þ 0 6 b¼6 4 fvz gx jti fpz gx jtþDt i 1 þ Dt qx Dz 3  2 7 7 fx qx fvz gx jitþDt 5 1 2qx DHx ðC 3Þ G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 C.2. Case II. Flow with a stationary bed 2 1 6 1 6 A ¼ 6 Dt 40 0 0 1 qx D z 0 0 3 fvz gx jtþDt i 6 fp g jtþDt 7 6 z xi 7 x¼6 7; 4 fvz gg jtþDt 5 i tþDt fpz gg ji fvz gx jti 1 6 6 fv g jt fp g jt 6 z xiþ z xi 6 b ¼ 6 Dt qx Dz 6 6 4 fvz gg jti 1 fpz gg jti 1 ðC 4Þ ðC 5Þ 0 1 2 2 3 0 0 0 07 7 7 1 05 1
2 fx fvz gx jti 1 2DHg Dzlb fvz gg jti Kbgzz 1 þ 3 2 7 7 7 7 2DHxg 7  2 7 7 t t fxg qx fvz gx ji 1 fvz gg ji 1 5  fxg fvz gx jti 1 fvz gg jti 1 2DHgx ðC 6Þ References Azar, J.J., Sanchez, R.A., 1997. Important Issues in Cuttings Transport for Drilling Directional Wells. Fifth Latin American and Caribbean Petroleum Engineering Conference and Exhibition, Rio de Janeiro, Brazil, Aug. 30–Sept. 3. SPE Paper 39020. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2002. Transport Phenomena, Second edition. John Wiley & Sons, Inc., New York. Carbonell, R.G., Whitaker, S., 1984. In: Bear, J., Corapcioglu, M.Y. (Eds.), Heat and Mass Transport in Porous Media. Mechanics of Fluid in Porous Media. Martinus Nijhoff, Bruselas. Cho, H., Shah, S.N., Osisanya, S.O., 2000. A three-layer modeling for cuttings transport with coiled tubing horizontal drilling. SPE Annual Technical Conference and Exhibition, Dallas, TX, October 1–4. SPE Paper 63269. Cho, H., Shah, S.N., Osisanya, S.O., 2002. A three-segment hydraulic model for cuttings transport in coiled tubing horizontal and deviated drilling. J. Can. Pet. Technol. 41 (6), 32–39. Doron, P., Barnea, D., 1993. A three-layer model for solid–liquid flow in horizontal pipes. Int. J. Multiph. Flow 19 (6), 1029–1043. Doron, P., Granica, D., Barnea, D., 1987. Slurry flow in horizontal pipes—experimental and modeling. Int. J. Multiph. Flow 13 (4), 535–547. Doron, P., Simkhis, M., Barnea, D., 1997. Flow of solid–liquid mixtures in inclined pipes. Int. J. Multiph. Flow 23 (2), 313–323. Drew, D.A., Cheng, L.Y., Lahey, R.T., 1979. The analysis of virtual mass effects in two phase flow. Int. J. Multiph. Flow 5, 233–242. Espinosa-Paredes, G., 2001. Theoretical derivation of the interaction effects whit an eccentric cell model and void fraction propagation in two-phase flow. Ann. Nucl. Energy 28, 659–688. Espinosa-Paredes, G., Cazarez-Candia, O., García-Gutiérrez, A., Martínez-Mendez, J., 2002. Void propagation in a bubbly twophase flow with expansion effects. Ann. Nucl. Energy 29, 1261–1298. 315 Ford, J.T., Peden, J.M., Oyeneyin, E.G., Zarrough, R., 1990. Experimental investigation of drilled cuttings transport in inclined boreholes. SPE Annual Technical Conference, New Orleans, Sept. 23–26. SPE Paper 20421. Gavignet, A.A., Sobey, I.J., 1989. Model aids cuttings transport prediction. JPT 916–921 (September). Gray, W.G., Lee, P.C.Y., 1977. On the theorems for local volume averaging of multiphase systems. Int. J. Multiph. Flow 3, 333–340. Huilin, L., Gidaspow, D., 2003. Hydrodynamics of binary fluidization in a riser: CFD simulation using two granular temperatures. Chem. Eng. Sci. 58, 3777–3792. Ishii, M., Mishima, K., 1984. Two-fluid model and hydrodynamic constitutive relations. Nucl. Eng. Des. 82, 107–126. Kamp, A.M., Rivero, M., 1999. Layer modeling for cuttings transport in highly inclined wellbores. SPE Latin American and Caribbean Petroleum Engineering Conference, Caracas, Venezuela, April 21– 23. SPE Paper 53942. Leising, L.J., Walton, I.C., 1998. Cuttings transport problems and solutions in coiled tubing drilling. IADC/SPE Drilling Conference, Dallas, TX, March 3–6. SPE Paper 39300. Li, J., Walker, S., 1999. Sensitive analysis of hole cleaning parameters in directional wells. SPE/ICoTA Coiled Tubing Roundtable, Houston, TX, May 25–26. SPE Paper 54498. Martins, A.L., Santana, C.C., 1992. Evaluation of cutting transport in horizontal and near horizontal wells—A dimensionless approach. SPE Latin American Petroleum Engineering Conference, Caracas, Venezuela, March 8–11. SPE Paper 23643. Martins, A.L., Sá, C.H.M., Lourenço, A.M.F., Freire, L.G.M., Campos, W., 1996. Experimental determination of interfacial friction factor in horizontal drilling with a bed of cuttings. SPE Latin American and Caribbean Petroleum Engineering Conference, Port of Spain, Trinidad and Tobago, April 23–26. SPE Paper 36075. Nguyen, D., Rahman, S.S., 1996. A three-layer hydraulic program for effective cuttings transport and hole cleaning in highly deviated and horizontal wells. IADC/SPE Asia Pacific Drilling Technology, Kuala Lumpur, Malaysia, Sept. 9–11. SPE Paper 36383. Nguyen, D., Rahman, S.S., 1998. A three-layer hydraulic program for effective cuttings transport and hole cleaning in highly deviated and horizontal wells. SPE Drilling and Completion, 182–189 (Sept.). SPE Paper 51186. Peden, J.M., Ford, J.T., Oyeneyin, M.B., 1990. Comprehensive experimental investigation of drilled cuttings transport in inclined wells including the effects of rotation and eccentricity. Europe 90, The Hague, Netherlands, October 22–24. SPE Paper 20925. Pilehvari, A.A., Azar, J.J., Shirazi, S.A., 1996. State-of the-art cuttings transport in horizontal wellbores. International Conference on Horizontal Well Technology, Calgary, Canada, Nov. 18–20. SPE Paper 37079. Salazar-Mendoza, M.R., 2004, Modelo hidrodinámico para el estudio del transporte de recortes de perforación en la sección horizontal de pozos petroleros, tesis de doctorado, Departamento de Mecánica, CENIDET, Cuernavaca, Mor., México. Sanchez, R.A., Azar, J.J., Bassal, A.A., Martins, A.L., 1999. Effect of drillpipe rotation on hole cleaning during directional well drilling. SPE J 4 (2), 101–107 (SPE Paper 56406). Santana, M., Martins, A.L., Sales, A., 1998. Advances in the modeling of the stratified flow of drilled cuttings in high angle and horizontal wells. SPE International Petroleum Conference and Exhibition, Villahermosa Tabasco, Mexico, March 3–5. SPE Paper 39890. 316 G. Espinosa-Paredes et al. / Journal of Petroleum Science and Engineering 55 (2007) 301–316 Sharma, M.P., 1990. Cuttings transport in inclined boreholes. Offshore South Asia Conference, Singapore, Dec. 4–7. OSEA Paper 90159. Sifferman, T.R., Becker, T.E., 1992. Hole cleaning in full-scale inclined wellbores. SPE Drilling Engineering, 115–120 (June). SPE Paper 20422. Televantos, Y., Shook, .A., Carleton, A., Streat, M., 1979. Flow of slurries of coarse particles at high solid concentrations. Can. J. Chem. Eng. 57, 255–262. Tomren, P.H., Iyoho, A.W., Azar, J.J., 1986. An experimental study of cuttings transport in directional wells. SPE Drilling Engineering, 43–56 (Feb.). SPE Paper 12123. Wallis, G.B., 1989. Inertial coupling in two-phase flow. Macroscopic properties of suspensions in an inviscid fluid. In: Hewitt, G.F., Delhaye, J.M., Zuber, N. (Eds.), Multiphase Science and Technology, vol 5. Hemisphere Publishing Corporation, New York, pp. 239–361. Whitaker, S., 1986. Flow in porous media I: A theoretical derivation of Darcy's Law. Transp. Porous Media 1, 3–25. Whitaker, S., 1996. The Forchheimer equation: a theoretical development. Transp. Porous Media 25, 27–61. Whitaker, S., 1999. The Method of Volume Averaging. Kluwer Academic Publishers, Dordrecht, The Netherlands. Wilson, K.C., 1987. Analysis of bed-load motion at high shear stress. J. Hydraul. Eng. 113 (1), 97–103 (ASCE Paper 21186). Zanotti, F., Carbonell, R.G., 1984. Development of transport equations for multiphase systems I: generalized development for two-phase systems. Chem. Eng. Sci. 39, 263–278.