- Email: [email protected]
Buoyancy driven slump flows of non-Newtonian fluids in pipes
Journal of Petroleum Science and Engineering 72 (2010) 236–243
Contents lists available at ScienceDirect
Journal of Petroleum Science and Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p e t r o l
Buoyancy driven slump flows of non-Newtonian fluids in pipes S. Malekmohammadi a, M.F. Naccache b, I.A. Frigaard a,c, D.M. Martinez d,⁎ a
Department of Mechanical Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z4 Department of Mechanical Engineering, Pontifícia Universidade Católica, Rio de Janeiro, RJ, Brazil 22453 c Department of Mathematics University of British Columbia, Vancouver, BC, Canada V6T 1Z2 d Department of Chemical and Biological Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z4 b
a r t i c l e
i n f o
Article history: Received 12 March 2009 Accepted 29 March 2010 Keywords: slump flows visco-plastic fluids yield stress exchange flows
a b s t r a c t We present results of a comparative experimental-numerical study of the buoyancy-driven slumping flow of two non-Newtonian fluids in a horizontal closed pipe. This type of flow occurs during the abandonment of horizontal oil and gas wells, when sealing the well through the setting of cement plugs. In our study, both of the fluids exhibit shear-thinning behaviour and one of the fluids has a yield stress. The yield stress is necessary in order for the flow to have a static equilibrium. We have studied the effects of changes in density difference and of small deviations from perfectly horizontal. The effects of these parameters on the interface shape and on the slump length versus time were analyzed. Comparison of numerical and experimental results shows broadly similar trends, but with some qualitative differences also observed, possibly due to interfacial effects at the interface. We also compare the slump length and shape against existing analytical results that predict the maximal static slump length. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved.
1. Introduction The motivation for this study comes from the industrial process of plug cementing, which is a common process in the oil industry; see e.g. Smith (1987); Nelson (1990). In this process dense cement slurry is placed into an oil or gas well that is already partially filled with another fluid. The aim of placing the cement plug is to form an impermeable hydraulic seal of the well and a hard mechanical barrier. A typical cement slurry density is around 1900 kg/m3, and the in-situ fluids, (drilling mud, viscous pill, spacer fluid, or simply residual reservoir fluids: oil/water), could have significantly lower density, say 1000– 1700 kg/m3. These density differences result in buoyancy induced fluid motion, which needs to be understood if the process is to be effective. Since the early 1990's there has been a massive increase in the numbers of oil wells constructed horizontally, primarily to increase productivity by aligning the well with the reservoir. The increase in horizontal wells was made possible by advances in directional drilling technology. The early 1990's saw a continual pushing of the horizontal reach of wells up to around 10 km, see e.g. the detailed description in Payne et al. (1995). The 10 km barrier was broken in a number of wells drilled at Wytch Farm, UK. The limits of “extreme” extended reach wells are now much higher, but these are not common. With
⁎ Corresponding author. E-mail address: [email protected] (D.M. Martinez).
present day technology it is feasible to reliably construct wells with horizontal extensions in the 7–10 km range. Whereas one of the initial drivers for horizontal wells was offshore multi-lateral drilling, reducing rig time and production footprint, these wells are now also common onshore. They occur both for conventional oil and gas wells and also in oil–sands development with the large-scale adoption of production processes such as SAGD (steam-assisted gravity drive) and others based on dissolving the heavy oil, typically with horizontal reaches of 1–3 km. Combined with the anticipated shift into increased gas production as world energy consumption climbs and oil reserves peak and decline over the next 20 years (Deffeyes, 2001; Goodstein, 2004; Tertzakian, 2006), issues of horizontal well abandonment will become increasingly topical. This is the focus of this paper. In abandoning a horizontal production well, we are typically placing a cement plug into low density Newtonian fluids. The buoyancy induced slumping motion can only be arrested by the yield stress of the cement slurry, which is typically insufficient to prevent initial motion during and after placement. Thus, slumping is inevitable and the yield stress acts to stop the flow only as gravitational forces are diminished by the slump extending down the pipe, reducing the slope of the interface. The final static state has been analysed by Frigaard and Ngwa (2004) who predicted both the maximum slump length and the shape of the slump. However, the slumping process is itself transient. In the initial transient phase of the flow inertial effects are potentially important and later there is the possibility of mixing, both of which may compromise the accuracy of the static slump estimates in Frigaard and Ngwa (2004). Thus, in this paper we investigate transient buoyancy-driven slumping flow inside a closed nearhorizontal pipe. We use both numerical and experimental methods,
0920-4105/$ – see front matter. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2010.03.024
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
as well as making comparison with the static predictions from Frigaard and Ngwa (2004). We present results showing the effects of density difference and pipe inclination on the slump length. We remark that there is a large literature on plug cementing in general and on this type of flow, e.g. (Beirute, 1978; Smith et al., 1983; Bour et al., 1986; Calvert et al., 1995; Harestad et al., 1997; Frigaard, 1998; Frigaard and Scherzer, 1998; Crawshaw and Frigaard, 1999; Fenie and Frigaard, 1999; Frigaard and Crawshaw, 1999; Frigaard and Scherzer, 2000). However, nearly all of these studies focus on plug-cementing in inclined-vertical wells, where the role of buoyancy is quite different. There have also been a large number of works in which thin-film/lubrication theory has been used to model the flow of a single visco-plastic fluid along a sloping surface, often inclined close to horizontal. In general the underlying momentum balance in these single fluid situations is much simpler to solve than in the 2-fluid case here. Consequently much progress has been made. For example, geophysical flows (lava, mudslides, etc) have been studied extensively; see Coussot (1997) and the reviews of Griffiths (2000); Balmforth et al. (2007). Other interesting studies of thin film flows of yield stress fluids are by Liu and Mei (1989); Piau (1996); Balmforth and Craster (1999); Ross et al. (2001); Wilson et al. (2002). An outline of our paper is as follows. Below in Section 2 we give an overview of the predictions of Frigaard and Ngwa (2004), on maximal static slump lengths and shapes, as it relates to our experiments and numerical modelling. Section 3 describes the methodologies that we have used in our study. In Section 4 we present our results, comparing first transient phenomena and then
237
against the static predictions. The paper concludes in Section 5 with a discussion of the principal results. 2. Analytical predictions In Frigaard and Ngwa (2004) the exchange flow of two Bingham fluids in a near-horizontal pipe, has been analyzed in the thin-film asymptotic limit. Although the transient problem is also formulated in Frigaard and Ngwa (2004), the principal results derived relate to the prediction of the maximal static slump length and shape. This maximal length is achieved when the shear stresses due to the density differences, slope of the interface and inclination of the pipe balance the yield stresses in each fluid. Since the fluid is static when the maximal length is achieved, the length of the slump depends only on the yield stresses and not on the effective viscosity of the fluids when yielded. Thus, the results in Frigaard and Ngwa (2004) are valid for any generalised Newtonian fluids with a yield stress. For the experimental situation we consider later, only one fluid (Carbopol) has a yield stress, τy, which simplifies the results. For our situation, the thin-film analysis in Frigaard and Ngwa (2004) requires that δ ≪ 1, where δ=
τy ; ½ρX −ρC gD
ð1Þ
D is the diameter of the pipe and ρX, ρC, denote respectively the densities of the Xanthan and Carbopol solutions. In our experiments
Fig. 1. Static analysis from Frigaard and Ngwa (2004): a) β0(h); b) βm(h); c) dimensionless slump shape for different χ χ. β0(h), βm(h), h and χ are all dimensionless parameters.
∊ [− 1, 1]; d) variation in dimensionless slump length with
238
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
we have δ typically between 0.03 and 0.1, so this requirement is approximately satisfied. The shape of the slump is described by the curve z(y), where y ∈ [0, D] denotes the height of the interface above the bottom of the pipe: dz D β0 ðhÞ ; = dh δ βm ðhÞ + χβ0 ðhÞ
ð2Þ
where h = y/D, where β0(h) and βm(h) are geometric functions defined in Frigaard and Ngwa (2004), (see Fig. 1a and b), and where χ = cosβ/δ. Here β denotes angle of well inclination from vertical, in degrees. The length of the slump, say Lstatic, is found from integrating the above: Lstatic =
D 1 β0 ðhÞ dh: ∫ δ 0 βm ðhÞ + χβ0 ðhÞ
ð3Þ
We note that dimensionless slump shapes and lengths are calculated by dividing through by D/δ = [ρX − ρC]gD2/τy. Fig. 1c plots the dimensionless slump shape for various χ ∈ [−1, 1]. For β b 90° we have positive χ. The slump is “uphill” and reduced by the inclination of the pipe. As we move to β N 90 we have χ N 0 and the effect of pipe inclination is to increase slump length. For the experiments conducted, values of χ lie within this interval. Fig. 1d plots the dimensionless slump length as a function of χ. 3. Methodology
Table 1 Series of fluids used in experiments. Series
1 2 3 4
Light fluid
Heavy fluid
Conc. Carbopol (wt.%)
Conc. Xanthan (wt.%)
Conc. sugar (wt.%)
Δρ/ρc
0.12 0.12 0.12 0.12
0.1 0.15 0.15 0.15
50 40 30 20
0.224 0.160 0.120 0.072
flat in the middle of the pipe. To minimize breakage of polymeric chains, pumps are driven with the lowest speed (5 Hz). The frame which holds the pipe is rotated to the horizontal position while the camera is capturing images continuously at a rate of 1 frame per second. Since most of slumping happens in the first 2 min, the camera speed is reduced gradually from 1 frame per second to 4 frames per hour. A mirror system allows views from the front and side of the pipe. The images are later processed to find the slump length in each timeframe. The length of slump is measured with a measuring tape that is mounted on the mirror system. The accuracy of measurements is 1 mm. To study the effect of pipe inclination 3 different angles were investigated: 89, 90 and 91°. Four sets of fluids were used in this work. Carbopol 940 was used in each set of experiment with constant aqueous concentration of 0.12% and a pH equal to 6.0. In order to see the effects of rheology and density difference 4 different Xanthan solutions were prepared, with
3.1. Experimental method The experimental set up consists of a long acrylic pipe with length of 1.83 m and inner diameter of D = 38 mm, an SLR digital camera, lighting system, two pumps, and a frame to support the pipe. A bearing connects the frame to the base steel structure and allows the pipe to be inclined from 0 to 93° (from vertical). Fig. 2 shows the schematic of the experimental set up. To perform the experiments the pipe is filled in the vertical position by a progressive cavity pump. The light fluid (Carbopol solution) is filled from the top and the heavy fluid (Xanthan solution), dyed black, is then filled from the bottom, leaving the initial interface
Fig. 2. Schematic of the experimental set-up.
Fig. 3. Typical flow curves: (a) 0.1% Xanthan solution at 23 °C. (b) 0.12% Carbopol 940 solution at pH 6.0 and 23 °C.
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
fluid properties, say ϕ, are shared by the fluids and represent volumeaveraged variables, defined by:
Table 2 Properties of the fluids used. Fluid 0.12% 0.10% 0.15% 0.15% 0.15%
Carbopol 940 Xanthan + 50% Xanthan + 40% Xanthan + 30% Xanthan + 20%
Sugar Sugar Sugar Sugar
239
Yield stress, τy (Pa)
Consistency, κ (Pa sn)
Power law index, n
Density ρ (kg/m3)
ϕ = α1 ϕ1 + α2 ϕ2
3–6 – – – –
21 0.6 0.3 0.2 0.2
0.6 0.7 0.6 0.7 0.5
1000 1224 1160 1120 1072
where ϕi is the solution variable or property of fluid i. Tracking of the fluids (hence the interface) is obtained by the solution of the continuity equation for phase 2:
ð4Þ
∂ðα2 ρ2 Þ + ∇⋅ðα2 ρ2 vÞ = 0 ∂t different concentrations of sugar. Sugar was used to increase both the density and viscosity of the Xanthan solutions. These fluids are described in Table 1. Δρ is the density difference between Heavy Fluid and Light Fluid and ρc is the density of Carbopol solution (Light Fluid). Note that of the rheological parameters, the yield stress is the most problematic to measure reliably and repeatably. By using the same Carbopol solution for each experiment and varying instead the Xanthan properties, we ensure consistency of any error, i.e. according to the static analysis the principal dimensionless parameter is the ratio of yield stress to buoyancy stress. The rheological properties of each solution were determined using a Bohlin C-VOR digital controlled shear stress-shear rate rheometer. Flow curves of a typical pair of experimental fluids, for 0.12% Carbopol 940 and 0.1% Xanthan + 50% sugar solutions, are shown in Fig. 3. Table 2 summarizes the measured rheological properties and density of each solution. The Levenberg–Marquardt method was employed for curve fitting of rheometric data. Xanthan solutions showed shear-thinning behaviour and were modeled by the powerlaw model. The Herschel–Bulkley model is commonly used to describe the rheological behavior of Carbopol solutions. However, determining the yield stress of a visco-plastic fluid is not an easy task. Fig. 4, suggests that there is a range of values for the yield stress rather than a single value. Although, there is not a single value for the yield stress, a constant value was needed to perform numerical simulation. As explained later, a value of 3.2 Pa was chosen for the yield stress of Carbopol in our numerical simulations. The rheological parameters of the fluids used are given in Table 2. 3.2. Numerical modeling The numerical solution was obtained using a finite volume method, (Patankar, 1980), with the 2 fluids modelled via the volume of fluid method (VOF) method. The VOF method solves a single set of momentum conservation equations in the domain and tracks the volume fraction of each fluid (αi, i = 1, 2). The solution variables and
Fig. 4. Range of yield stresses of Carbopol, using the differentiation method.
ð5Þ
where t is time, ρ2 is the density of fluid 2 and v is the velocity vector. Note that Eq. (5) also represents a concentration equation for the case in which molecular diffusion between fluids is vanishingly small. The volume fraction of fluid 1 is obtained by the condition: α1 + α2 = 1. The momentum equation is given by: h i ∂ðρvÞ T + ∇⋅ðρvvÞ = −∇p + ∇⋅ η ∇v + ∇v + ρg ∂t
ð6Þ
where ρ is the density, p is the pressure, g is the gravitational acceleration vector and η is the viscosity function. Again note that ρ and η need be interpolated from the properties of the individual phases. Thus, we see that the two fluids are handled via the well-known volume of fluid method, which is equivalent to assuming the two phases are in a mixture. At each point in space and time, a single velocity field, pressure field and volume fraction is computed and tracked: the interface is interpreted as the position where the volume fraction is equal to 0.5. For our numerical simulations we always denote the heavier fluid (Xanthan) by fluid 1. The rheological model for fluid 1 is the powerlaw model: n −1 η1 = κ1γ˙ 1
ð7Þ
(See Ostwald (1925) and Blair et al. (1939)). For fluid 2, (Carbopol), a modified bi-viscosity function was used to model the visco-plastic behavior: 8τ < y + κ γ˙ n2 −1 2 η2 = γ˙ : η0
9 if γ˙ N γ˙ small = otherwise
;
ð8Þ
where τy is the yield stress, η0 is a low shear rate viscosity constant and γ̇small = τy/η0. Here η1, κ1 and n1 are viscosity, consistency and power law index for fluid 1, respectively. η2, κ2, n2 are viscosity, consistency and power law index for fluid 2, respectively. Note that since the units of η1 and η2 are Pa s and n is a constant, κ1 and κ2 have units of Pa sn. Apart from an initial stage of the flow, inertial effects are minimal and although the flow is 3D the assumption of symmetry about a vertical plane cutting the pipe along the centerline of the pipe is reasonable. Boundary conditions are the usual no-slip and impermeability conditions at the solid boundaries, and symmetry conditions on the central plane of the pipe. In the experiments the tube was filled with both fluids in a vertical position and then rotated to the horizontal position. To mimic this, at t= 0 the (numerical) tube was in vertical position with interface a horizontal interface which was perpendicular to the direction of gravity. For tN 0 the numerical domain was rotated to horizontal at constant angular speed during the initial 3 s of the simulation. Although inertial effects were included in all simulations, after the initial few seconds the inertial terms in the equations are negligible. Numerical results (Fig. 5) showed that inertia effects were only important in the beginning of the simulations and did not affect the eventual evolution of the slump length. The tube diameter in the computations was identical to that in the experiments, D = 38 mm, but the length was about 1/3 that in the experiments, (i.e. 0.6 m), in order to reduce computational times. Computational tests were run to ensure that the end walls of the tube
240
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
Fig. 5. Inertial effects in a typical simulation.
were not affecting flow at the interface. Mesh tests were performed on 3 unstructured meshes: Meshes 1–3 with 124,080, 60,900 and 26,474 control volumes, respectively. The differences in slump length over time between Mesh 2 and 3 were up to 10%; those between Mesh 1 and 2 were below 5%. Three-dimensional computations with complex fluids are costly and those on Mesh 1 took approximately 10 times longer than those on Mesh 2. Thus as a compromise between speed and accuracy, Mesh 2 was used for our comparative results. The software Fluent (Hirt and Nichols, 1981; Fluent User's Guide, 2006) is used for all computations.
4. Results 4.1. Transient results The effects of pipe inclination and density difference between fluids on the slump length were investigated. The rheological parameters used in the numerical simulations for the lighter fluid (Carbopol solution) were: η0 = 10,000 Pa s, τy = 3.2 Pa, n = 0.7, and κ = 21 Pa sn. To choose a reasonable value for yield stress within the range indicated in Fig. 4, different values for the yield stress were tried in the numerical simulations. The numerical results of different simulations suggested the value of τy = 3.2 Pa gives the closest results to the experimental results for all 4 series. With reference to Fig. 4, we observe that this value is also close to the maximum. As mentioned before, four different heavier fluids (Xanthan solutions) were analyzed, each modeled by the Power-law equation. The rheological parameters used in the simulations were as those given in Table 2.
Fig. 7. (a) Experimental results for the slump length versus time for series 1, Δρ/ρc = 0.224. (b) Numerical results for the slump length versus time for series 1, Δρ/ρc = 0.224.
Fig. 6 shows the interface shapes after 1.5 h, in series 1 and series 2 of experiments, at 90° inclination of pipe (Heavier fluid is the dark one). It can be observed that for higher density differences the slump length is longer and the interface shape is more stretched, as expected. The flow visualization shows an irregular interface, which could suggest some possible mixing close to the interface. A reasonable agreement between numerical and experimental results is observed.
Fig. 6. Comparison of interface shape after 1.5 h: (a) Experimental result: series 1, Δρ/ρc =0.224; (b) Numerical result: series 1, Δρ/ρc =0.224; (c) Experimental result: series 2, Δρ/ρc =0.160; (d) Numerical result: series 2, Δρ/ρc =0.160.
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
241
Fig. 7 shows (a) experimental and (b) numerical results of the dimensionless slump length increase with time, for three pipe inclinations. The dimensionless slump length and time are defined as:
T
L =
L D
t T t = h 1 = n i n+1 D = ΔρgD =κ
ð9Þ
The experimental errors for a typical case are shown in Fig. 8. The highest and lowest values of slump length are shown with dots in this figure. The solid line is the curve fitting of average values of slump lengths. It is observed that there is a steep increase in the slump length at the beginning of the process, but as flow goes on the slump length tends to an asymptote. At the onset, the buoyancy driven stresses are high enough to overcome the yield stress, but as the interface stretches buoyancy effects decrease and eventually the flow stops. We also observe a slight increase in slump length as the pipe inclination varies from 89° to 91°. Some differences between numerical and experimental results are observed. The numerical solution gives lower slump lengths, but a higher slump velocity at later times. This behavior could be explained by a mixing process that may be going on, changing the actual fluids properties close to the interface, or could be due to the modeling methodology: either the use of a viscosity regularization technique or the multiphase modeling using VOF. Moreover, it is noted that the experimental errors are larger for lower times, where the comparison is worse. 4.2. Quasi-static results Fig. 9a and b shows the slump length as a function of the density difference for three pipe inclinations, after 1.5 h. Fig. 9c shows the slump length predicted by the analytical method. A linear relation is observed, for the analytical, numerical and experimental results. It is noted that the qualitative behavior between these results is similar, but some discrepancy in the quantitative values is observed. Since the slump flow is buoyancy driven, it is to be expected that pipe inclination effects become more important at higher density differences. 5. Discussion The results we have presented show reasonable quantitative and qualitative agreement between numerical and experimental results. However, there are still obvious discrepancies that need explaining.
Fig. 9. (a) Experimental results for the slump length after 1.5 h versus Δρ/ρc for α = 89°, 90° and 91°; (b) Numerical results for the slump length after 1.5 h versus Δρ/ρc for α = 89°, 90° and 91°; (c) Analytical results for the slump length versus Δρ/ρc for α = 89°, 90° and 91°.
Fig. 8. Experimental results for the slump length versus time for series 1, Δρ/ρc = 0.224, α = 90°.
Apart from purely numerical artifacts, we feel that two different physical explanations could lead to the discrepancies observed between the numerical and experimental results. First note that,
242
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
after an initial phase in which inertia may be significant, we expect that the chief physical balance in the flow is between buoyancy and yield stress, (i.e. viscous stresses are decaying as the flow decelerates, as are inertial stresses). Thus, the first potential source of error comes from the low shear rheology. Measurement of low shear rate rheology is a difficult task and, in the case of yield stress fluids, there is a ongoing debate about the true nature of the material behaviour at sub-yield stress values. In essence the question is whether the sub-yield behaviour of a given fluid is better described by slow viscous flow or by elastic deformation. For the purpose of numerical models such as Fluent, the former of these behaviours is assumed. Low shear rheological effects are governed by both τy and η0: as the shear stress drops below τy the fluid behaves as a very viscous fluid with viscosity η0. The effect of low shear rate viscosity on slump length through time is shown in Fig. 10. It is observed that by increasing η0 ten times, slump length does not change initially. However, for larger times, the velocities obtained for very low η0 fluid are much larger than the ones observed experimentally, where the flow became quasi static after almost 2 h. It is also observed that for η0 N 5000 Pa s, the results seem to be unaffected by η0. Thus, taking η0 = 10,000 Pa s is a reasonable approximation for η0. As mentioned in Section 4.1, in our numerical simulations we have used a value of 3.2 Pa close to the lower limit of the range of experimental yield stress values (3 Pa). This is partly because of the possibility of mixing. If little mixing happens very close to the interface, it would cause a decrease of the effective yield stress of the lighter fluid, and also a decrease of the effective density difference between the two fluids. Measuring the rheological properties and densities of the fluids at the interface during and after the experiments is not possible. However the effect of this assumption can be investigated in the numerical simulations. Fig. 11 shows how decreasing the yield stress of the lighter fluid affects the slump length. It is observed that with a very high value, τy = 6.0 Pa, the discrepancy between experimental and numerical results is large. However, as we decrease the yield stress, the numerical results approach the experimental ones. It is worth mentioning that at lower times the experimental errors are larger. The second possible explanation for the discrepancy is related to possible mixing of the two fluids close to the interface, during the experiments (see Fig. 6). This mixing could both reduce the local yield stress, as discussed above, and lower the density difference between fluids. There is notably less mixing in the simulations than the experiment, essentially due to the numerical VOF model which
Fig. 10. Experimental and numerical results for the slump length versus time, for series 1, Δρ/ρc = 0.224 and α = 90°. Effect of the lighter fluid's low shear rate viscosity.
Fig. 11. Experimental and numerical results for the slump length versus time, for series 1, Δρ/ρc = 0.224 and α = 90°. Effect of the lighter fluid's yield stress.
attempts to keep the interface sharp. We may however simulate this mixing effect by artificially decreasing the density difference between the fluids after some fixed (arbitrary) time in numerical simulation. In doing so, lower velocities were obtained at larger times and a similar final transient behavior between numerical and experimental results could be observed; see Fig. 12. Being able to get a very close fit between experimental and numerical results this way, does support the idea of mixing at the interface being responsible, but the fitting method is evidently somewhat ad hoc. Possibly this could be improved by introducing a proper miscible fluids formulation. It should however be mentioned that mixing cannot progress quickly. Carbopol solution is a gel and consists of large molecular networks which hinder the mixing of the fluids. In the experimental results shown earlier we have seen an apparently asymptotic increase in the slump length towards some finite limit. It is worth asking if there is some kind of simplified model in which the same characteristics are found. One option here is to attempt to derive and solve a thin-film style model, for which the results in Section 2 provide the static limit. This is made difficult by the pipe geometry, but would be feasible for a plane channel geometry. A transient model of this type has been examined in Fenie and Frigaard (1999), for planar channels that are inclined close to the vertical. In
Fig. 12. Experimental and numerical results for the slump length versus time, for series 1 and α = 90°. Effect of the variation of Δρ/ρc in time.
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
that model it is possible to show that there is a finite stopping time. There are in fact a large number of examples in the literature of such finite time decay results, e.g. (Glowinski, 1983; Huilgol et al., 2002). For strongly inclined channels, the effects of the interface slope are neglected at leading order. This leads to a hyperbolic system for the interface motion. We are thus able to bound the rate of spreading of the interface and derive estimates of this kind. For a strictly horizontal pipe (or other duct), there is no driving force for the fluid motion apart from the buoyancy gradient generated through the slope of the interface. This results in a parabolic problem for the interface height and the velocity is locally a function of the interface slope as well as its height. Recently, in Matson and Hogg (2007, 2008) it has been shown that for thin-film motions driven by the interface slope there is no finite stopping time. Instead, in approaching the final static steady state the flow variables decay like t− n, with n denoting the power law index. The flow considered in Matson and Hogg (2007, 2008) is the classical thin film flow, for a Herschel–Bulkley fluid. Here we have a more complex geometry and a 2-fluid flow satisfying a zero flux condition. However, the driving force for the flow is still buoyancy rather than a forced flux and we believe that a similar result should hold. For both our experimental and numerical results, we do not appear to have a finite time stopping. We did not collect extensive data at long times, but approximately the slump length appears to asymptote to a constant value, with discrepancy decaying algebraically, like t− ϕ where ϕ ∈ (1, 2), but varies between experiments. In terms of practical consequences of our results, perhaps the main one is to confirm that the static slump length estimates computed in Frigaard and Ngwa (2004), are reasonable, even if conservative. Fig. 9 shows very similar trends between these estimates, numerical simulations and experimental lengths, over a range of hole inclinations. The main method for using such estimates would be to calculate a surplus volume of cement, to take care of the slump volume. This volume is essentially wastage since no hydraulic seal is achieved in the slump. Secondly, our study has demonstrated that computational fluid dynamics (CFD) tools, such as those utilised, can be effective in giving a close approximation to lab experiments. Where this becomes important is in situations where the transient phase is important. Although it is the yield stress versus buoyancy balance that dictates the final length of the slump, it is the effective viscosity which controls the transient timescale. CFD is thus a useful tool for experimenting with different fluid rheologies, in viscous pills, spacers and the like. Evidently, there are large benefits in terms of cost, compared with the running of lab experiments or yard tests. Finally, we need to discuss some of the limitations of our study. In the oil industry a case study approach is becoming increasingly common in studying technical problems areas. For a process such as the setting of abandonment plugs, such an approach has the advantage of being able to include the complexities of the real process, e.g. cement setting, temperature effects, salinity of reservoir fluids, effects of placement methodology, etc. These complexities are those that we have eliminated from our study in order to get a grasp on the fluid dynamics of post placement slumping. Although neglecting these effects does introduces limitations, we note that the initial fluid motions (responsible for a large part of the slump length) occur fast, over a timescale of seconds/minutes, whereas most of the effects neglected do take place on significantly longer timescales. This explains our reasons for neglecting them. Some disadvantages of the case study approach over a controlled/ scientific approach used here are the variability of methods used for setting plugs, the difficulty in getting access to good quality data from the wellsite, and particularly in the case of abandonment plugs, the
243
lack of rigorous evaluation of whether or not plugs have actually sealed the well.
Acknowledgments This research was supported by NSERC Canada (SM, IAF & DMM) and byCNPq Brazil (MFN).
References Balmforth, N.J., Craster, R.V., 1999. A consistent thin-layer theory for Bingham plastics. J. Non-Newton. Fluid Mech. 84, 65–81. Balmforth, N.J., Craster, R.V., Rust, A.C., Sassi, R., 2007. Viscoplastic flow over an inclined surface. J. Non-Newton. Fluid Mech. 142, 219–243. Beirute, R.M., 1978. Flow behaviour of an Unset Cement Plug in Place. SPE paper 7589. Blair, G.W.S., Hening, J.C., Wagstaff, A., 1939. The flow of cream through narrow glass tubes. J. Phys. Chem. 43 (7), 853–864. Bour, D.L., Sutton, D.L., Creel, P.G., 1986. Development of Effective Methods for Placing Competent Cement Plugs. SPE paper 15008. Calvert, D.G., Heathman, J.F., Griffith, J.E., 1995. Plug Cementing: Horizontal to Vertical Conditions. SPE paper 30514. Coussot, P., 1997. Mud flow Rheology and Dynamics. IAHR/AIRH Monograph. Crawshaw, J.P., Frigaard, I., 1999. Cement Plugs: Stability and Failure by BuoyancyDriven Mechanism. SPE paper 56959. Deffeyes, K.S., 2001. Hubbert's peak : the Impending World Oil Shortage. Princeton University Press. Fenie, H., Frigaard, I.A., 1999. Preventing buoyancy driven flows of two Bingham fluids in a closed pipe: fluid rheology design for oilfield plug cementing. Math. Comput. Model. 30 (7–8), 71–91. Fluent User's Guide, 2006. v. 6.2. Fluent Inc., New Hampshire. Frigaard, I.A., 1998. Stratified exchange flows of two Bingham fluids in an inclined slot. J. Non-Newton. Fluid Mech. 78, 61–87. Frigaard, I.A., Crawshaw, J.P., 1999. Preventing buoyancy driven flows of two Bingham fluids in a closed pipe: Fluid rheology design for oilfield plug. J. Eng. Math. 36 (4), 327–348. Frigaard, I.A., Ngwa, N.G., 2004. Upper bounds on the slump length in plug cementing of near-horizontal wells. J. Non-Newton. Fluid Mech. 117 (2–3), 147–162. Frigaard, I.A., Scherzer, O., 1998. Uniaxial exchange flows of two Bingham fluids in a cylindrical duct. IMA J. Appl. Math. 61, 237–266. Frigaard, I.A., Scherzer, O., 2000. The effects of yield stress variation on uniaxial exchange flows of two Bingham fluids in a pipe. SIAM J. Appl. Math. 60 (6), 1950–1976. Glowinski, R., 1983. Numerical Methods for Nonlinear Variational Problems. SpringerVerlag. Goodstein, D.L., 2004. Out of gas : the end of the age of oil. W.W. Norton & Company Inc. Griffiths, R.W., 2000. The dynamics of lava flows. Annu. Rev. Fluid Mech. 32, 477–518. Harestad, K., Herigstad, T.P., Torsvoll, A., Nodland, N.E., Saasen, A., 1997. Optimization of Balanced-Plug Cementing. SPE paper 33084. Hirt, C.W., Nichols, B.D., 1981. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 204–225. Huilgol, R.R., Mena, B., Piau, J.M., 2002. Finite stopping time problems and rheometry of Bingham fluids. J. Non-Newton. Fluid Mech. 102, 97–107. Liu, K.F., Mei, C.C., 1989. Slow spreading of a sheet of Bingham fluid on an inclined plane. J. Fluid Mech. 207, 505–529. Matson, G.P., Hogg, A.J., 2007. Two-dimensional dam break flows of Herschel–Bulkley fluids: the approach to the arrested state. J. Non-Newton. Fluid Mech. 142, 79–94. Matson, G.P., Hogg, A.J., 2008. Slumps of viscoplastic fluids on slopes. submitted to J. Non-Newton. Fluid Mech. Nelson, E.B., 1990. Well Cementing. Schlumberger Educational Services. Ostwald, W., 1925. The velocity function of viscosity of disperse systems. Kolloid Z. 36, 99–117. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation. Payne, M.L., Wilton, B.S., Ramos, G.G., 1995. Recent Advances and Emerging Technologies for Extended Reach Drilling. SPE paper 29920. Piau, J.M., 1996. Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane. J. Rheol. 40, 711–723. Ross, A.B., Wilson, S.K., Duffy, B.R., 2001. Thin-film flow of a viscoplastic material round a large horizontal stationary or rotating cylinder. J. Fluid Mech. 430, 309–333. Smith, D.K., 1987. Cementing. SPE monograph series, Society of Petroleum Engineers. Smith, R.C., Beirute, R.M., Holman, G.B., 1983. Improved Method of Setting Successful Whipstock Cement Plugs. SPE paper 11415. Tertzakian, P., 2006. A thousand barrels a second : the coming oil break point and the challenges facing an energy dependent world. McGraw-Hill. Wilson, S.K., Duffy, B.R., Ross, A.B., 2002. On the gravity-driven draining of a rivulet of viscoplastic material down a slowly varying substrate. Phys. Fluids 14, 555–571.
Contents lists available at ScienceDirect
Journal of Petroleum Science and Engineering j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / p e t r o l
Buoyancy driven slump flows of non-Newtonian fluids in pipes S. Malekmohammadi a, M.F. Naccache b, I.A. Frigaard a,c, D.M. Martinez d,⁎ a
Department of Mechanical Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z4 Department of Mechanical Engineering, Pontifícia Universidade Católica, Rio de Janeiro, RJ, Brazil 22453 c Department of Mathematics University of British Columbia, Vancouver, BC, Canada V6T 1Z2 d Department of Chemical and Biological Engineering, University of British Columbia, Vancouver, BC, Canada V6T 1Z4 b
a r t i c l e
i n f o
Article history: Received 12 March 2009 Accepted 29 March 2010 Keywords: slump flows visco-plastic fluids yield stress exchange flows
a b s t r a c t We present results of a comparative experimental-numerical study of the buoyancy-driven slumping flow of two non-Newtonian fluids in a horizontal closed pipe. This type of flow occurs during the abandonment of horizontal oil and gas wells, when sealing the well through the setting of cement plugs. In our study, both of the fluids exhibit shear-thinning behaviour and one of the fluids has a yield stress. The yield stress is necessary in order for the flow to have a static equilibrium. We have studied the effects of changes in density difference and of small deviations from perfectly horizontal. The effects of these parameters on the interface shape and on the slump length versus time were analyzed. Comparison of numerical and experimental results shows broadly similar trends, but with some qualitative differences also observed, possibly due to interfacial effects at the interface. We also compare the slump length and shape against existing analytical results that predict the maximal static slump length. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved.
1. Introduction The motivation for this study comes from the industrial process of plug cementing, which is a common process in the oil industry; see e.g. Smith (1987); Nelson (1990). In this process dense cement slurry is placed into an oil or gas well that is already partially filled with another fluid. The aim of placing the cement plug is to form an impermeable hydraulic seal of the well and a hard mechanical barrier. A typical cement slurry density is around 1900 kg/m3, and the in-situ fluids, (drilling mud, viscous pill, spacer fluid, or simply residual reservoir fluids: oil/water), could have significantly lower density, say 1000– 1700 kg/m3. These density differences result in buoyancy induced fluid motion, which needs to be understood if the process is to be effective. Since the early 1990's there has been a massive increase in the numbers of oil wells constructed horizontally, primarily to increase productivity by aligning the well with the reservoir. The increase in horizontal wells was made possible by advances in directional drilling technology. The early 1990's saw a continual pushing of the horizontal reach of wells up to around 10 km, see e.g. the detailed description in Payne et al. (1995). The 10 km barrier was broken in a number of wells drilled at Wytch Farm, UK. The limits of “extreme” extended reach wells are now much higher, but these are not common. With
⁎ Corresponding author. E-mail address: [email protected] (D.M. Martinez).
present day technology it is feasible to reliably construct wells with horizontal extensions in the 7–10 km range. Whereas one of the initial drivers for horizontal wells was offshore multi-lateral drilling, reducing rig time and production footprint, these wells are now also common onshore. They occur both for conventional oil and gas wells and also in oil–sands development with the large-scale adoption of production processes such as SAGD (steam-assisted gravity drive) and others based on dissolving the heavy oil, typically with horizontal reaches of 1–3 km. Combined with the anticipated shift into increased gas production as world energy consumption climbs and oil reserves peak and decline over the next 20 years (Deffeyes, 2001; Goodstein, 2004; Tertzakian, 2006), issues of horizontal well abandonment will become increasingly topical. This is the focus of this paper. In abandoning a horizontal production well, we are typically placing a cement plug into low density Newtonian fluids. The buoyancy induced slumping motion can only be arrested by the yield stress of the cement slurry, which is typically insufficient to prevent initial motion during and after placement. Thus, slumping is inevitable and the yield stress acts to stop the flow only as gravitational forces are diminished by the slump extending down the pipe, reducing the slope of the interface. The final static state has been analysed by Frigaard and Ngwa (2004) who predicted both the maximum slump length and the shape of the slump. However, the slumping process is itself transient. In the initial transient phase of the flow inertial effects are potentially important and later there is the possibility of mixing, both of which may compromise the accuracy of the static slump estimates in Frigaard and Ngwa (2004). Thus, in this paper we investigate transient buoyancy-driven slumping flow inside a closed nearhorizontal pipe. We use both numerical and experimental methods,
0920-4105/$ – see front matter. Crown Copyright © 2010 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2010.03.024
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
as well as making comparison with the static predictions from Frigaard and Ngwa (2004). We present results showing the effects of density difference and pipe inclination on the slump length. We remark that there is a large literature on plug cementing in general and on this type of flow, e.g. (Beirute, 1978; Smith et al., 1983; Bour et al., 1986; Calvert et al., 1995; Harestad et al., 1997; Frigaard, 1998; Frigaard and Scherzer, 1998; Crawshaw and Frigaard, 1999; Fenie and Frigaard, 1999; Frigaard and Crawshaw, 1999; Frigaard and Scherzer, 2000). However, nearly all of these studies focus on plug-cementing in inclined-vertical wells, where the role of buoyancy is quite different. There have also been a large number of works in which thin-film/lubrication theory has been used to model the flow of a single visco-plastic fluid along a sloping surface, often inclined close to horizontal. In general the underlying momentum balance in these single fluid situations is much simpler to solve than in the 2-fluid case here. Consequently much progress has been made. For example, geophysical flows (lava, mudslides, etc) have been studied extensively; see Coussot (1997) and the reviews of Griffiths (2000); Balmforth et al. (2007). Other interesting studies of thin film flows of yield stress fluids are by Liu and Mei (1989); Piau (1996); Balmforth and Craster (1999); Ross et al. (2001); Wilson et al. (2002). An outline of our paper is as follows. Below in Section 2 we give an overview of the predictions of Frigaard and Ngwa (2004), on maximal static slump lengths and shapes, as it relates to our experiments and numerical modelling. Section 3 describes the methodologies that we have used in our study. In Section 4 we present our results, comparing first transient phenomena and then
237
against the static predictions. The paper concludes in Section 5 with a discussion of the principal results. 2. Analytical predictions In Frigaard and Ngwa (2004) the exchange flow of two Bingham fluids in a near-horizontal pipe, has been analyzed in the thin-film asymptotic limit. Although the transient problem is also formulated in Frigaard and Ngwa (2004), the principal results derived relate to the prediction of the maximal static slump length and shape. This maximal length is achieved when the shear stresses due to the density differences, slope of the interface and inclination of the pipe balance the yield stresses in each fluid. Since the fluid is static when the maximal length is achieved, the length of the slump depends only on the yield stresses and not on the effective viscosity of the fluids when yielded. Thus, the results in Frigaard and Ngwa (2004) are valid for any generalised Newtonian fluids with a yield stress. For the experimental situation we consider later, only one fluid (Carbopol) has a yield stress, τy, which simplifies the results. For our situation, the thin-film analysis in Frigaard and Ngwa (2004) requires that δ ≪ 1, where δ=
τy ; ½ρX −ρC gD
ð1Þ
D is the diameter of the pipe and ρX, ρC, denote respectively the densities of the Xanthan and Carbopol solutions. In our experiments
Fig. 1. Static analysis from Frigaard and Ngwa (2004): a) β0(h); b) βm(h); c) dimensionless slump shape for different χ χ. β0(h), βm(h), h and χ are all dimensionless parameters.
∊ [− 1, 1]; d) variation in dimensionless slump length with
238
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
we have δ typically between 0.03 and 0.1, so this requirement is approximately satisfied. The shape of the slump is described by the curve z(y), where y ∈ [0, D] denotes the height of the interface above the bottom of the pipe: dz D β0 ðhÞ ; = dh δ βm ðhÞ + χβ0 ðhÞ
ð2Þ
where h = y/D, where β0(h) and βm(h) are geometric functions defined in Frigaard and Ngwa (2004), (see Fig. 1a and b), and where χ = cosβ/δ. Here β denotes angle of well inclination from vertical, in degrees. The length of the slump, say Lstatic, is found from integrating the above: Lstatic =
D 1 β0 ðhÞ dh: ∫ δ 0 βm ðhÞ + χβ0 ðhÞ
ð3Þ
We note that dimensionless slump shapes and lengths are calculated by dividing through by D/δ = [ρX − ρC]gD2/τy. Fig. 1c plots the dimensionless slump shape for various χ ∈ [−1, 1]. For β b 90° we have positive χ. The slump is “uphill” and reduced by the inclination of the pipe. As we move to β N 90 we have χ N 0 and the effect of pipe inclination is to increase slump length. For the experiments conducted, values of χ lie within this interval. Fig. 1d plots the dimensionless slump length as a function of χ. 3. Methodology
Table 1 Series of fluids used in experiments. Series
1 2 3 4
Light fluid
Heavy fluid
Conc. Carbopol (wt.%)
Conc. Xanthan (wt.%)
Conc. sugar (wt.%)
Δρ/ρc
0.12 0.12 0.12 0.12
0.1 0.15 0.15 0.15
50 40 30 20
0.224 0.160 0.120 0.072
flat in the middle of the pipe. To minimize breakage of polymeric chains, pumps are driven with the lowest speed (5 Hz). The frame which holds the pipe is rotated to the horizontal position while the camera is capturing images continuously at a rate of 1 frame per second. Since most of slumping happens in the first 2 min, the camera speed is reduced gradually from 1 frame per second to 4 frames per hour. A mirror system allows views from the front and side of the pipe. The images are later processed to find the slump length in each timeframe. The length of slump is measured with a measuring tape that is mounted on the mirror system. The accuracy of measurements is 1 mm. To study the effect of pipe inclination 3 different angles were investigated: 89, 90 and 91°. Four sets of fluids were used in this work. Carbopol 940 was used in each set of experiment with constant aqueous concentration of 0.12% and a pH equal to 6.0. In order to see the effects of rheology and density difference 4 different Xanthan solutions were prepared, with
3.1. Experimental method The experimental set up consists of a long acrylic pipe with length of 1.83 m and inner diameter of D = 38 mm, an SLR digital camera, lighting system, two pumps, and a frame to support the pipe. A bearing connects the frame to the base steel structure and allows the pipe to be inclined from 0 to 93° (from vertical). Fig. 2 shows the schematic of the experimental set up. To perform the experiments the pipe is filled in the vertical position by a progressive cavity pump. The light fluid (Carbopol solution) is filled from the top and the heavy fluid (Xanthan solution), dyed black, is then filled from the bottom, leaving the initial interface
Fig. 2. Schematic of the experimental set-up.
Fig. 3. Typical flow curves: (a) 0.1% Xanthan solution at 23 °C. (b) 0.12% Carbopol 940 solution at pH 6.0 and 23 °C.
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
fluid properties, say ϕ, are shared by the fluids and represent volumeaveraged variables, defined by:
Table 2 Properties of the fluids used. Fluid 0.12% 0.10% 0.15% 0.15% 0.15%
Carbopol 940 Xanthan + 50% Xanthan + 40% Xanthan + 30% Xanthan + 20%
Sugar Sugar Sugar Sugar
239
Yield stress, τy (Pa)
Consistency, κ (Pa sn)
Power law index, n
Density ρ (kg/m3)
ϕ = α1 ϕ1 + α2 ϕ2
3–6 – – – –
21 0.6 0.3 0.2 0.2
0.6 0.7 0.6 0.7 0.5
1000 1224 1160 1120 1072
where ϕi is the solution variable or property of fluid i. Tracking of the fluids (hence the interface) is obtained by the solution of the continuity equation for phase 2:
ð4Þ
∂ðα2 ρ2 Þ + ∇⋅ðα2 ρ2 vÞ = 0 ∂t different concentrations of sugar. Sugar was used to increase both the density and viscosity of the Xanthan solutions. These fluids are described in Table 1. Δρ is the density difference between Heavy Fluid and Light Fluid and ρc is the density of Carbopol solution (Light Fluid). Note that of the rheological parameters, the yield stress is the most problematic to measure reliably and repeatably. By using the same Carbopol solution for each experiment and varying instead the Xanthan properties, we ensure consistency of any error, i.e. according to the static analysis the principal dimensionless parameter is the ratio of yield stress to buoyancy stress. The rheological properties of each solution were determined using a Bohlin C-VOR digital controlled shear stress-shear rate rheometer. Flow curves of a typical pair of experimental fluids, for 0.12% Carbopol 940 and 0.1% Xanthan + 50% sugar solutions, are shown in Fig. 3. Table 2 summarizes the measured rheological properties and density of each solution. The Levenberg–Marquardt method was employed for curve fitting of rheometric data. Xanthan solutions showed shear-thinning behaviour and were modeled by the powerlaw model. The Herschel–Bulkley model is commonly used to describe the rheological behavior of Carbopol solutions. However, determining the yield stress of a visco-plastic fluid is not an easy task. Fig. 4, suggests that there is a range of values for the yield stress rather than a single value. Although, there is not a single value for the yield stress, a constant value was needed to perform numerical simulation. As explained later, a value of 3.2 Pa was chosen for the yield stress of Carbopol in our numerical simulations. The rheological parameters of the fluids used are given in Table 2. 3.2. Numerical modeling The numerical solution was obtained using a finite volume method, (Patankar, 1980), with the 2 fluids modelled via the volume of fluid method (VOF) method. The VOF method solves a single set of momentum conservation equations in the domain and tracks the volume fraction of each fluid (αi, i = 1, 2). The solution variables and
Fig. 4. Range of yield stresses of Carbopol, using the differentiation method.
ð5Þ
where t is time, ρ2 is the density of fluid 2 and v is the velocity vector. Note that Eq. (5) also represents a concentration equation for the case in which molecular diffusion between fluids is vanishingly small. The volume fraction of fluid 1 is obtained by the condition: α1 + α2 = 1. The momentum equation is given by: h i ∂ðρvÞ T + ∇⋅ðρvvÞ = −∇p + ∇⋅ η ∇v + ∇v + ρg ∂t
ð6Þ
where ρ is the density, p is the pressure, g is the gravitational acceleration vector and η is the viscosity function. Again note that ρ and η need be interpolated from the properties of the individual phases. Thus, we see that the two fluids are handled via the well-known volume of fluid method, which is equivalent to assuming the two phases are in a mixture. At each point in space and time, a single velocity field, pressure field and volume fraction is computed and tracked: the interface is interpreted as the position where the volume fraction is equal to 0.5. For our numerical simulations we always denote the heavier fluid (Xanthan) by fluid 1. The rheological model for fluid 1 is the powerlaw model: n −1 η1 = κ1γ˙ 1
ð7Þ
(See Ostwald (1925) and Blair et al. (1939)). For fluid 2, (Carbopol), a modified bi-viscosity function was used to model the visco-plastic behavior: 8τ < y + κ γ˙ n2 −1 2 η2 = γ˙ : η0
9 if γ˙ N γ˙ small = otherwise
;
ð8Þ
where τy is the yield stress, η0 is a low shear rate viscosity constant and γ̇small = τy/η0. Here η1, κ1 and n1 are viscosity, consistency and power law index for fluid 1, respectively. η2, κ2, n2 are viscosity, consistency and power law index for fluid 2, respectively. Note that since the units of η1 and η2 are Pa s and n is a constant, κ1 and κ2 have units of Pa sn. Apart from an initial stage of the flow, inertial effects are minimal and although the flow is 3D the assumption of symmetry about a vertical plane cutting the pipe along the centerline of the pipe is reasonable. Boundary conditions are the usual no-slip and impermeability conditions at the solid boundaries, and symmetry conditions on the central plane of the pipe. In the experiments the tube was filled with both fluids in a vertical position and then rotated to the horizontal position. To mimic this, at t= 0 the (numerical) tube was in vertical position with interface a horizontal interface which was perpendicular to the direction of gravity. For tN 0 the numerical domain was rotated to horizontal at constant angular speed during the initial 3 s of the simulation. Although inertial effects were included in all simulations, after the initial few seconds the inertial terms in the equations are negligible. Numerical results (Fig. 5) showed that inertia effects were only important in the beginning of the simulations and did not affect the eventual evolution of the slump length. The tube diameter in the computations was identical to that in the experiments, D = 38 mm, but the length was about 1/3 that in the experiments, (i.e. 0.6 m), in order to reduce computational times. Computational tests were run to ensure that the end walls of the tube
240
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
Fig. 5. Inertial effects in a typical simulation.
were not affecting flow at the interface. Mesh tests were performed on 3 unstructured meshes: Meshes 1–3 with 124,080, 60,900 and 26,474 control volumes, respectively. The differences in slump length over time between Mesh 2 and 3 were up to 10%; those between Mesh 1 and 2 were below 5%. Three-dimensional computations with complex fluids are costly and those on Mesh 1 took approximately 10 times longer than those on Mesh 2. Thus as a compromise between speed and accuracy, Mesh 2 was used for our comparative results. The software Fluent (Hirt and Nichols, 1981; Fluent User's Guide, 2006) is used for all computations.
4. Results 4.1. Transient results The effects of pipe inclination and density difference between fluids on the slump length were investigated. The rheological parameters used in the numerical simulations for the lighter fluid (Carbopol solution) were: η0 = 10,000 Pa s, τy = 3.2 Pa, n = 0.7, and κ = 21 Pa sn. To choose a reasonable value for yield stress within the range indicated in Fig. 4, different values for the yield stress were tried in the numerical simulations. The numerical results of different simulations suggested the value of τy = 3.2 Pa gives the closest results to the experimental results for all 4 series. With reference to Fig. 4, we observe that this value is also close to the maximum. As mentioned before, four different heavier fluids (Xanthan solutions) were analyzed, each modeled by the Power-law equation. The rheological parameters used in the simulations were as those given in Table 2.
Fig. 7. (a) Experimental results for the slump length versus time for series 1, Δρ/ρc = 0.224. (b) Numerical results for the slump length versus time for series 1, Δρ/ρc = 0.224.
Fig. 6 shows the interface shapes after 1.5 h, in series 1 and series 2 of experiments, at 90° inclination of pipe (Heavier fluid is the dark one). It can be observed that for higher density differences the slump length is longer and the interface shape is more stretched, as expected. The flow visualization shows an irregular interface, which could suggest some possible mixing close to the interface. A reasonable agreement between numerical and experimental results is observed.
Fig. 6. Comparison of interface shape after 1.5 h: (a) Experimental result: series 1, Δρ/ρc =0.224; (b) Numerical result: series 1, Δρ/ρc =0.224; (c) Experimental result: series 2, Δρ/ρc =0.160; (d) Numerical result: series 2, Δρ/ρc =0.160.
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
241
Fig. 7 shows (a) experimental and (b) numerical results of the dimensionless slump length increase with time, for three pipe inclinations. The dimensionless slump length and time are defined as:
T
L =
L D
t T t = h 1 = n i n+1 D = ΔρgD =κ
ð9Þ
The experimental errors for a typical case are shown in Fig. 8. The highest and lowest values of slump length are shown with dots in this figure. The solid line is the curve fitting of average values of slump lengths. It is observed that there is a steep increase in the slump length at the beginning of the process, but as flow goes on the slump length tends to an asymptote. At the onset, the buoyancy driven stresses are high enough to overcome the yield stress, but as the interface stretches buoyancy effects decrease and eventually the flow stops. We also observe a slight increase in slump length as the pipe inclination varies from 89° to 91°. Some differences between numerical and experimental results are observed. The numerical solution gives lower slump lengths, but a higher slump velocity at later times. This behavior could be explained by a mixing process that may be going on, changing the actual fluids properties close to the interface, or could be due to the modeling methodology: either the use of a viscosity regularization technique or the multiphase modeling using VOF. Moreover, it is noted that the experimental errors are larger for lower times, where the comparison is worse. 4.2. Quasi-static results Fig. 9a and b shows the slump length as a function of the density difference for three pipe inclinations, after 1.5 h. Fig. 9c shows the slump length predicted by the analytical method. A linear relation is observed, for the analytical, numerical and experimental results. It is noted that the qualitative behavior between these results is similar, but some discrepancy in the quantitative values is observed. Since the slump flow is buoyancy driven, it is to be expected that pipe inclination effects become more important at higher density differences. 5. Discussion The results we have presented show reasonable quantitative and qualitative agreement between numerical and experimental results. However, there are still obvious discrepancies that need explaining.
Fig. 9. (a) Experimental results for the slump length after 1.5 h versus Δρ/ρc for α = 89°, 90° and 91°; (b) Numerical results for the slump length after 1.5 h versus Δρ/ρc for α = 89°, 90° and 91°; (c) Analytical results for the slump length versus Δρ/ρc for α = 89°, 90° and 91°.
Fig. 8. Experimental results for the slump length versus time for series 1, Δρ/ρc = 0.224, α = 90°.
Apart from purely numerical artifacts, we feel that two different physical explanations could lead to the discrepancies observed between the numerical and experimental results. First note that,
242
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
after an initial phase in which inertia may be significant, we expect that the chief physical balance in the flow is between buoyancy and yield stress, (i.e. viscous stresses are decaying as the flow decelerates, as are inertial stresses). Thus, the first potential source of error comes from the low shear rheology. Measurement of low shear rate rheology is a difficult task and, in the case of yield stress fluids, there is a ongoing debate about the true nature of the material behaviour at sub-yield stress values. In essence the question is whether the sub-yield behaviour of a given fluid is better described by slow viscous flow or by elastic deformation. For the purpose of numerical models such as Fluent, the former of these behaviours is assumed. Low shear rheological effects are governed by both τy and η0: as the shear stress drops below τy the fluid behaves as a very viscous fluid with viscosity η0. The effect of low shear rate viscosity on slump length through time is shown in Fig. 10. It is observed that by increasing η0 ten times, slump length does not change initially. However, for larger times, the velocities obtained for very low η0 fluid are much larger than the ones observed experimentally, where the flow became quasi static after almost 2 h. It is also observed that for η0 N 5000 Pa s, the results seem to be unaffected by η0. Thus, taking η0 = 10,000 Pa s is a reasonable approximation for η0. As mentioned in Section 4.1, in our numerical simulations we have used a value of 3.2 Pa close to the lower limit of the range of experimental yield stress values (3 Pa). This is partly because of the possibility of mixing. If little mixing happens very close to the interface, it would cause a decrease of the effective yield stress of the lighter fluid, and also a decrease of the effective density difference between the two fluids. Measuring the rheological properties and densities of the fluids at the interface during and after the experiments is not possible. However the effect of this assumption can be investigated in the numerical simulations. Fig. 11 shows how decreasing the yield stress of the lighter fluid affects the slump length. It is observed that with a very high value, τy = 6.0 Pa, the discrepancy between experimental and numerical results is large. However, as we decrease the yield stress, the numerical results approach the experimental ones. It is worth mentioning that at lower times the experimental errors are larger. The second possible explanation for the discrepancy is related to possible mixing of the two fluids close to the interface, during the experiments (see Fig. 6). This mixing could both reduce the local yield stress, as discussed above, and lower the density difference between fluids. There is notably less mixing in the simulations than the experiment, essentially due to the numerical VOF model which
Fig. 10. Experimental and numerical results for the slump length versus time, for series 1, Δρ/ρc = 0.224 and α = 90°. Effect of the lighter fluid's low shear rate viscosity.
Fig. 11. Experimental and numerical results for the slump length versus time, for series 1, Δρ/ρc = 0.224 and α = 90°. Effect of the lighter fluid's yield stress.
attempts to keep the interface sharp. We may however simulate this mixing effect by artificially decreasing the density difference between the fluids after some fixed (arbitrary) time in numerical simulation. In doing so, lower velocities were obtained at larger times and a similar final transient behavior between numerical and experimental results could be observed; see Fig. 12. Being able to get a very close fit between experimental and numerical results this way, does support the idea of mixing at the interface being responsible, but the fitting method is evidently somewhat ad hoc. Possibly this could be improved by introducing a proper miscible fluids formulation. It should however be mentioned that mixing cannot progress quickly. Carbopol solution is a gel and consists of large molecular networks which hinder the mixing of the fluids. In the experimental results shown earlier we have seen an apparently asymptotic increase in the slump length towards some finite limit. It is worth asking if there is some kind of simplified model in which the same characteristics are found. One option here is to attempt to derive and solve a thin-film style model, for which the results in Section 2 provide the static limit. This is made difficult by the pipe geometry, but would be feasible for a plane channel geometry. A transient model of this type has been examined in Fenie and Frigaard (1999), for planar channels that are inclined close to the vertical. In
Fig. 12. Experimental and numerical results for the slump length versus time, for series 1 and α = 90°. Effect of the variation of Δρ/ρc in time.
S. Malekmohammadi et al. / Journal of Petroleum Science and Engineering 72 (2010) 236–243
that model it is possible to show that there is a finite stopping time. There are in fact a large number of examples in the literature of such finite time decay results, e.g. (Glowinski, 1983; Huilgol et al., 2002). For strongly inclined channels, the effects of the interface slope are neglected at leading order. This leads to a hyperbolic system for the interface motion. We are thus able to bound the rate of spreading of the interface and derive estimates of this kind. For a strictly horizontal pipe (or other duct), there is no driving force for the fluid motion apart from the buoyancy gradient generated through the slope of the interface. This results in a parabolic problem for the interface height and the velocity is locally a function of the interface slope as well as its height. Recently, in Matson and Hogg (2007, 2008) it has been shown that for thin-film motions driven by the interface slope there is no finite stopping time. Instead, in approaching the final static steady state the flow variables decay like t− n, with n denoting the power law index. The flow considered in Matson and Hogg (2007, 2008) is the classical thin film flow, for a Herschel–Bulkley fluid. Here we have a more complex geometry and a 2-fluid flow satisfying a zero flux condition. However, the driving force for the flow is still buoyancy rather than a forced flux and we believe that a similar result should hold. For both our experimental and numerical results, we do not appear to have a finite time stopping. We did not collect extensive data at long times, but approximately the slump length appears to asymptote to a constant value, with discrepancy decaying algebraically, like t− ϕ where ϕ ∈ (1, 2), but varies between experiments. In terms of practical consequences of our results, perhaps the main one is to confirm that the static slump length estimates computed in Frigaard and Ngwa (2004), are reasonable, even if conservative. Fig. 9 shows very similar trends between these estimates, numerical simulations and experimental lengths, over a range of hole inclinations. The main method for using such estimates would be to calculate a surplus volume of cement, to take care of the slump volume. This volume is essentially wastage since no hydraulic seal is achieved in the slump. Secondly, our study has demonstrated that computational fluid dynamics (CFD) tools, such as those utilised, can be effective in giving a close approximation to lab experiments. Where this becomes important is in situations where the transient phase is important. Although it is the yield stress versus buoyancy balance that dictates the final length of the slump, it is the effective viscosity which controls the transient timescale. CFD is thus a useful tool for experimenting with different fluid rheologies, in viscous pills, spacers and the like. Evidently, there are large benefits in terms of cost, compared with the running of lab experiments or yard tests. Finally, we need to discuss some of the limitations of our study. In the oil industry a case study approach is becoming increasingly common in studying technical problems areas. For a process such as the setting of abandonment plugs, such an approach has the advantage of being able to include the complexities of the real process, e.g. cement setting, temperature effects, salinity of reservoir fluids, effects of placement methodology, etc. These complexities are those that we have eliminated from our study in order to get a grasp on the fluid dynamics of post placement slumping. Although neglecting these effects does introduces limitations, we note that the initial fluid motions (responsible for a large part of the slump length) occur fast, over a timescale of seconds/minutes, whereas most of the effects neglected do take place on significantly longer timescales. This explains our reasons for neglecting them. Some disadvantages of the case study approach over a controlled/ scientific approach used here are the variability of methods used for setting plugs, the difficulty in getting access to good quality data from the wellsite, and particularly in the case of abandonment plugs, the
243
lack of rigorous evaluation of whether or not plugs have actually sealed the well.
Acknowledgments This research was supported by NSERC Canada (SM, IAF & DMM) and byCNPq Brazil (MFN).
References Balmforth, N.J., Craster, R.V., 1999. A consistent thin-layer theory for Bingham plastics. J. Non-Newton. Fluid Mech. 84, 65–81. Balmforth, N.J., Craster, R.V., Rust, A.C., Sassi, R., 2007. Viscoplastic flow over an inclined surface. J. Non-Newton. Fluid Mech. 142, 219–243. Beirute, R.M., 1978. Flow behaviour of an Unset Cement Plug in Place. SPE paper 7589. Blair, G.W.S., Hening, J.C., Wagstaff, A., 1939. The flow of cream through narrow glass tubes. J. Phys. Chem. 43 (7), 853–864. Bour, D.L., Sutton, D.L., Creel, P.G., 1986. Development of Effective Methods for Placing Competent Cement Plugs. SPE paper 15008. Calvert, D.G., Heathman, J.F., Griffith, J.E., 1995. Plug Cementing: Horizontal to Vertical Conditions. SPE paper 30514. Coussot, P., 1997. Mud flow Rheology and Dynamics. IAHR/AIRH Monograph. Crawshaw, J.P., Frigaard, I., 1999. Cement Plugs: Stability and Failure by BuoyancyDriven Mechanism. SPE paper 56959. Deffeyes, K.S., 2001. Hubbert's peak : the Impending World Oil Shortage. Princeton University Press. Fenie, H., Frigaard, I.A., 1999. Preventing buoyancy driven flows of two Bingham fluids in a closed pipe: fluid rheology design for oilfield plug cementing. Math. Comput. Model. 30 (7–8), 71–91. Fluent User's Guide, 2006. v. 6.2. Fluent Inc., New Hampshire. Frigaard, I.A., 1998. Stratified exchange flows of two Bingham fluids in an inclined slot. J. Non-Newton. Fluid Mech. 78, 61–87. Frigaard, I.A., Crawshaw, J.P., 1999. Preventing buoyancy driven flows of two Bingham fluids in a closed pipe: Fluid rheology design for oilfield plug. J. Eng. Math. 36 (4), 327–348. Frigaard, I.A., Ngwa, N.G., 2004. Upper bounds on the slump length in plug cementing of near-horizontal wells. J. Non-Newton. Fluid Mech. 117 (2–3), 147–162. Frigaard, I.A., Scherzer, O., 1998. Uniaxial exchange flows of two Bingham fluids in a cylindrical duct. IMA J. Appl. Math. 61, 237–266. Frigaard, I.A., Scherzer, O., 2000. The effects of yield stress variation on uniaxial exchange flows of two Bingham fluids in a pipe. SIAM J. Appl. Math. 60 (6), 1950–1976. Glowinski, R., 1983. Numerical Methods for Nonlinear Variational Problems. SpringerVerlag. Goodstein, D.L., 2004. Out of gas : the end of the age of oil. W.W. Norton & Company Inc. Griffiths, R.W., 2000. The dynamics of lava flows. Annu. Rev. Fluid Mech. 32, 477–518. Harestad, K., Herigstad, T.P., Torsvoll, A., Nodland, N.E., Saasen, A., 1997. Optimization of Balanced-Plug Cementing. SPE paper 33084. Hirt, C.W., Nichols, B.D., 1981. Volume of Fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 204–225. Huilgol, R.R., Mena, B., Piau, J.M., 2002. Finite stopping time problems and rheometry of Bingham fluids. J. Non-Newton. Fluid Mech. 102, 97–107. Liu, K.F., Mei, C.C., 1989. Slow spreading of a sheet of Bingham fluid on an inclined plane. J. Fluid Mech. 207, 505–529. Matson, G.P., Hogg, A.J., 2007. Two-dimensional dam break flows of Herschel–Bulkley fluids: the approach to the arrested state. J. Non-Newton. Fluid Mech. 142, 79–94. Matson, G.P., Hogg, A.J., 2008. Slumps of viscoplastic fluids on slopes. submitted to J. Non-Newton. Fluid Mech. Nelson, E.B., 1990. Well Cementing. Schlumberger Educational Services. Ostwald, W., 1925. The velocity function of viscosity of disperse systems. Kolloid Z. 36, 99–117. Patankar, S.V., 1980. Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation. Payne, M.L., Wilton, B.S., Ramos, G.G., 1995. Recent Advances and Emerging Technologies for Extended Reach Drilling. SPE paper 29920. Piau, J.M., 1996. Flow of a yield stress fluid in a long domain. Application to flow on an inclined plane. J. Rheol. 40, 711–723. Ross, A.B., Wilson, S.K., Duffy, B.R., 2001. Thin-film flow of a viscoplastic material round a large horizontal stationary or rotating cylinder. J. Fluid Mech. 430, 309–333. Smith, D.K., 1987. Cementing. SPE monograph series, Society of Petroleum Engineers. Smith, R.C., Beirute, R.M., Holman, G.B., 1983. Improved Method of Setting Successful Whipstock Cement Plugs. SPE paper 11415. Tertzakian, P., 2006. A thousand barrels a second : the coming oil break point and the challenges facing an energy dependent world. McGraw-Hill. Wilson, S.K., Duffy, B.R., Ross, A.B., 2002. On the gravity-driven draining of a rivulet of viscoplastic material down a slowly varying substrate. Phys. Fluids 14, 555–571.