Concentric cylinder viscometer flows of Herschel-Bulkley fluids

Abstract: Drilling fluids and well cements are example non-Newtonian fluids that are used for geothermal and petroleum well construction. Measurement of the nonNewtonian fluid viscosities are normally performed using a concentric cylinder Couette geometry, where one of the cylinders rotates at a controlled speed or under a controlled torque. In this paper we address Couette flow of yield stress shear thinningfluids in concentric cylinder geometries.We focus on typical oilfield viscometers and discuss effects of yield stress and shear thinning on fluid yielding at low viscometer rotational speeds and errors caused by the Newtonian shear rate assumption. We relate these errors to possible implications for typical wellbore flows.


Introduction
In the enterprise of constructing geothermal or petroleum wells, drilling fluids and well cements are used. These fluids are constructed as a blend of particle suspensions, liquid-liquid dispersions and polymer solutions with a large variety of additional chemicals. The fluids show a strong non-Newtonian characteristic.
As an example, in Figure 1 we plot typical shapes of flow curves for fluids involved in the operation where a casing string is cemented to the newly drilled formation. Primary cementing is performed once the current well section is drilled to target depth and a casing string is run into the hole and to the bottom of the section. The operation schedule involves drilling, where the drilling fluid transports out drilled cuttings, displacement of the drilling fluid from the narrow annular space between the casing string and the formation and replace it by a cement slurry to create a total zonal isolation in the well. Once in place, the cement slurry is allowed to harden into a cement sheath. A successful operation is dependent on the complete displacement of drilling fluid from the narrow annular space outside the casing. This is typically achieved by injecting a sequence of wash and spacer fluids ahead of the cement slurry, and the displacement is optimized through careful design of fluid densities and viscosities. As illustrated by the example flow curves in Figure 1, the involved fluids typically exhibit yield stress and thus show shear thinning behavior. Most drilling fluids, and some cements slurries also exhibit a curved shape of the viscosity flow curve after reaching the yield stress; i.e. yield stress shear thinning fluids as proposed by Herschel and Bulkley [1]. In addition to exhibiting shear thinning and yield stress behaviour, most drilling and well construction fluids are also thixotropic, i.e. requiring a certain time duration for microstructure equilibration to current shearing conditions [2,3]. In this paper, we neglect thixotropic behaviour and instead focus on the steady state, non-Newtonian flow curve of well construction fluids. Stable displacements are normally preferred for primary cementing operations in the laminar regime, implying that the displacing fluid express a higher frictional pressure loss than the drilling fluid to optimise displacement. In practice this means that the cement should be effectively more viscous than the fluid it displaces [4,5]. It is also essential to limit the circulating pressures in the well to avoid hydraulic fracturing the formation and uncontrolled loss of fluids. This means that to optimize displacement and avoid formation fracturing during pumping, it is important to reduce the viscosity of the drilling fluids being in the annulus. However, this reduction cannot be so large that the pressure control in the well is lost. Therefore, it is important to have accurate and reliable viscosity characterization of the fluids involved in the processes.
In the field, fluid viscosities are normally measured using a concentric cylinder Couette geometry, where the outer cup rotates at a fixed angular velocity Ω and the inner bob is held fixed by application of a torque M. The oilfield standard for viscosity measurements are based on equipment description. The dimensions for the conventional oilfield viscometer were described by Savins and Roper [6]. The rotor-and-bob geometry is as follows: Rotor radius is 1.8415 cm and the bob radius is 1.7247 cm. The length of the bob is 3.8 cm. These dimensions were partly sought to present the viscosity at the shear rate of 511 s −1 directly in units of centipoise and to allow for measurements of fluids containing weighting particles of sizes less than 75 micron. According to oilfield standards for well cementing, the viscometer should have at least the following rotation speeds: 600, 300, 200, 100, 6 and 3 revolutions per minute (RPM). The 600 RPM is used for drilling fluids and spacers, but not for well cements as exposure to this shear rate reduces the reproducibility of flow curves. Later, it has become an increasing use of viscometers with more rotation speeds in the oilfields. Typically, additional rotation speeds equal to 60, 30, 20, 10, 2, and 1 RPM can now occasionally be found on oilfield equipment. The dimensions remain the same.
In the absence of end effects, the applied torque is related to the wall shear stress at the bob via where the submerged bob radius and length are denoted R 1 and L, respectively. This results in a relationship between angular velocity of the cup and the corresponding shear stress at the bob surface. However, a relationship between the shear rate and the corresponding shear stress is typically what is sought for.
For a concentric cylinder viscometer, Krieger and Elrod [7] showed that the angular velocity Ω of the rotating cup is related to the shear rate and shear stress in the gap via Eq. (2) where τ 1 is again the wall shear stress at the inner bob and κ = R 1 /R 2 is the ratio of bob to cup radii. As Ω is controlled and τ 1 is measured, Eq. (2) is to be solved for the unknown shear rate˙(τ 1 ). Krieger and Elrod [7] and later Krieger [8] showed that Eq. (2) can be solved for the unknown shear rate via an infinite series without making a priori assumptions of the type of flow curve for the fluid. These solutions do not always perform well for fluids with yield stresses [9] as yield stress introduces a mathematical singularity in the viscosity curve. Alternative solution strategies based on Tikhonov regularization [10,11] and wavelet-vaguelette decomposition [12] have been developed, which are capable of handling yield stress fluids with no a priori flow curve assumption.
The problem of existence of yield stresses is summarised by Watson [13]. However, for our purpose we consider sufficiently short time frames where use of an yield stress is applicable. Therefore, as will be described later, we will assume an existence of the yield stresses.
If the relationship between shear rate and shear stress is known a priori for the fluid in question, relations linking angular velocity, applied torque and model parameters can be derived from Eq. (2), see e.g. Li et al. [14]. Following a different approach, Chatzimina et al. [15] solved the governing equations for Herschel-Bulkley fluids and found expressions for the wall shear rate in a concentric cylinder geometry. They showed that wall shear rates may deviate significantly from corresponding Newtonian values depending on the diameter ratio of the Couette geometry and material parameters.
In this paper, we assume a priori yield stress shear thinning (Herschel-Bulkley) flow curves and study the effects of yield stress and shear thinning on the wall shear rates. We compare the actual Herschel-Bulkley wall shear rates to those of a Newtonian fluid and quantify the error in the common Newtonian assumption for model yield stress fluids. The outline of the paper is as follows: We begin by considering the momentum equation for the Herschel-Bulkley fluid in a Couette geometry in the next section. We derive an equation relating rotational velocity of the cup to the shear stress in the gap that is equivalent to what can be obtained through Eq. (2). Next we solve the governing equation for azimuthal flow velocity in the gap and bob shear rate for a certain, fixed Herschel-Bulkley fluid. Finally, we consider example Fann 35 model viscometer measurements of well construction fluids, including a typical water-based drilling fluid, a spacer fluid and cement slurry, and investigate the difference between the Newtonian shear rate assumption and the actual non-Newtonian wall shear rate. We briefly discuss implications before concluding.

Couette flow equations
We consider a concentric cylinder Couette geometry where the outer cylinder rotates at a constant angular velocity Ω and the inner bob is held fixed by application of a constant torque. We focus on Herschel-Bulkley fluids described by the constitutive equation: Since the length of the viscometer gap is very much longer than the width of the gap, we neglect end-effects and assume the fluid velocity in the gap is purely azimuthal and only a function of the radial position in the gap, ⃗ v = v θ (r)⃗ e θ . Subject to this assumption, the shear rate is given by˙= and the azimuthal momentum equation is or τ rθ = τ(r) = c/r 2 where c is a constant of integration. We follow the derivation presented by Chatzimina et al. [15] and use thaṫ = r ∂ ∂r in the sheared region in the annular gap between the two cylinders. If the outer cylinder rotates sufficiently fast, or if the yield stress is sufficiently small, the fluid will be fully sheared throughout the gap between the cylinders. Denote the minimum rotational velocity of the outer cup that ensures fully sheared fluid by Ω * . This value can be obtained by choosing c such that τ(R 2 ) = τy or equivalently˙= 0 at r = R 2 . Integrating Eq. (6) with the boundary condition v θ (R 2 ) = ΩR 2 one finds: We evaluate the integral in Eq. (7) numerically for different aspect ratios κ and three different shear thinning indices in Figure 2. As expected, the function decreases rapidly as κ increases and the gap between the bob and the cup narrows.
For Ω ≥ Ω * , the fluid is considered fully sheared in the gap between the bob and the cup and the wall shear rate at the bob is given by Eq. (6) with r = R 1 and c determined implicitly from Eq. (8) (︂ kΩ n τy The left hand side of Eq. (8) is the ratio of viscous stress to yield stress in the fluid to the power 1/n. This ratio is often expressed as the Bingham number, i.e. Bn = τy /(kΩ n ), as per Chatzimina et al. [15]. In the case where τy → 0, n → 1 and k = µ, we regain the Newtonian result where c = 2µΩκ 2 R 2 2 /(1 − κ 2 ) and the wall shear rate at the bob is˙w ,N = 2Ω/(1 − κ 2 ), i.e. only dependent on rotational velocity Ω and the aspect ratio of the cup, κ. For a shear thinning power law fluid, τ = k˙n, the wall shear rate can also be solved analytically, with the result being˙w ,PL = (2Ω/n)/(1 − κ (2/n) ) [16]. As there is an analytical relation between wall shear rate and rotational velocity Ω, it is possible to fit measurements of shear stress as function of Ω to the constitutive model via τ = k[(2Ω/n)/(1 − κ (2/n) )] n without invoking a Newtonian assumption about the shear rate.
The existence of a non-zero yield stress for Herschel-Bulkley fluids makes the analysis more complicated than the Newtonian and power law results above, since there now may be a region of unyielded fluid in the gap; the existence of unyielded fluid is dependent on the model pa-rameters to be determined. If Ω < Ω * , or equivalently if the Bingham number is larger than a critical Bingham number, the fluid is sheared only in the gap between R 1 and R plug ; a rigid plug region exists between R plug and the outer cylinder wall at R 2 . In this case, τ(r) = τy R 2 plug /r 2 , and R plug is determined implicitly from (︂ kΩ n τy with κp = R plug /R 2 . The wall shear rate at the bob is then obtained from Eq. (4), as before. partly sheared fluid and a plug region next to the outer cup that moves as a rigid body. The other speeds shown in Figure 3 result in fully sheared fluid. As seen from the slope of the normalized angular velocity at the bob surface, the wall shear rates for the Herschel-Bulkley fluid exceed that of a Newtonian fluid at all rotational speeds. We will compare wall shear rates in more detail below.

Wall shear rates
In Figures 4 to 6 we plot the ratio of actual wall shear rates to the commonly assumed Newtonian shear rate at the inner bob wall as function of rotational speed of the outer cup for different combinations of the Herschel-Bulkley parameters. The Fann Model 35 speeds are marked on the curves by crosses. In all cases, the difference between the actual shear rate and the Newtonian assumption increases with increasing degree of non-Newtonian behaviour; that is, decreasing n, increasing τy and decreas-  ing k enhance differences between the two wall shear rates. The wall shear rate ratio is a monotonically decreasing function of rotational speed of the cup; it is at the lowest rotational speeds that non-Newtonian flow curve most affects the velocity profile and correspondingly the wall shear rate at the bob.

Example application
As specific examples, we consider Fann Model 35 R1-B1 measurements of different well fluids; a water-based drilling fluid (WBM), a water-based spacer and a cement slurry. The density of the drilling fluid is 1250 kg/m 3 , the spacer density is 1700 kg/m 3 and the cement slurry density is 1920 kg/m 3 . Fann Model 35 measurements of the three samples are listed in Table 1. We note that the measurements listed in the table correspond to measured angular deflections at different rotational speeds of the outer cup. The measured angles are converted to wall shear stress by multiplying the angle measurement by the convesion factor 0.5107 Pa/degree. The measurements are acquired using the R1-B1 measurement geometry which consists of an inner bob with radius 0.017245 m and an outer cup of radius 0.018415 m. The ratio of inner to outer radius is then κ ≈ 0.94.
We now proceed to fit these measurements to the Herschel-Bulkley model using the results presented in the previous sections in order to obtain the actual wall shear rate at the bob. We perform this numerically by searching for the Herschel-Bulkley model parameter combination that minimizes the sum of squared residuals using the Levenberg-Marquardt method [17]. Every time the model   Table 1. We plot fits assuming both standard Newtonian shear rates, power law shear rates and Herschel-Bulkley shear rates.
parameters are changed during the iterative search, we first evaluate Eq. (7) to obtain the critical rotational speed below which the fluid is partly yielded. We secondly estimate the actual bob wall shear rates corresponding to each of the six rotational speeds using Eq. (6) with the integration constant in τ determined using either Eq. (8) or (9).

Water-based drilling fluid; WBM
We begin our analysis by considering the WBM measurements in Table 1. In Figure 7 we plot as points the Fann Model 35 measurements assuming either Newtonian, power law or Herschel-Bulkley wall shear rates at the bob. As anticipated for a shear thinning fluid, the actual wall shear rates are larger than the assumed Newtonian shear rates. The solid lines correspond to the best Herschel-Bulkley model fits for the two data sets obtained using the method outlined above. Using the Newtonian shear rate assumption, the best approximation to the data points is found to be the flow curve τ N = (4.56 + 0.73˙0 .52 ) Pa.
For power law shear rates, we obtain the parametrization τ PL = (4.56+0.70˙0 .52 ) Pa, while when utilizing the actual Herschel-Bulkley wall shear rates, we find the flow curve τ HB = (4.29 + 0.70˙0 .52 ) Pa. In other words, both Newtonian and power law wall shear rates results in a small overestimation of the fluid yield stress compared to when using Herschel-Bulkley wall shear rates. The power law approximation however results in very similar numerical values for k and n as for Herschel-Bulkley shear rates.
In Table 2 we list the Newtonian shear rates in a Fann Model 35 R1-B1 measurement geometry as well as the estimated power law and Herschel-Bulkley shear rates for the above parametrization. As discussed above, the Newtonian shear rate is given by the rotational speed of the cup and the aspect ratio,˙w ,N = 2Ω/(1 − κ 2 ) while the power law shear rate at the bob is˙w ,PL = (2Ω/n)/(1 − κ (2/n) ). The Herschel-Bulkley shear rate is obtained from Eq. (4).
As one would expect, the power law wall shear rate approximation produces shear rate estimates intermediate between the Newtonian and the Herschel-Bulkley estimates. The power law approximation approaches that of the Herschel-Bulkley result as the rotational speed of the cup increases, as expected. The power law shear rate approximation simplifies the numerical processing compared to the Herschel-Bulkley solution presented above, as the power law wall shear rate is given explicitly through the aspect ratio, Ω and n. The Herschel-Bulkley solution is implicit, requiring more computations to evaluate the existence of plug region and determination of the shear stress in the gap, Eqs. (8) and (9).

Spacer fluid
Next, we analyse the spacer measurements in Table 1. Using the Newtonian shear rate assumption on this data set, the best approximation to the data points is found to be  Table 1. We plot fits assuming both standard Newtonian shear rates, power law shear rates and Herschel-Bulkley shear rates. the flow curve τ N = (6.16 + 0.80˙0 .54 ) Pa. For power law shear rates, we obtain the parametrization τ PL = (6.16 + 0.78˙0 .54 ) Pa, while when utilizing the actual Herschel-Bulkley wall shear rates, we find the flow curve τ HB = (5.79 + 0.78˙0 .54 ) Pa. As per the WBM in the previous section, both Newtonian and power law wall shear rates result in small over-estimations of the yield stress compared to the Herschel-Bulkley wall shear rates. As before, the power law approximation however results in very similar numerical values for k and n as for Herschel-Bulkley shear rates. We plot the Fann 35 measurements and the curve fits for the spacer fluid in Figure 8. For completeness, the Newtonian, power law and Herschel-Bulkley shear rates in a Fann Model 35 R1-B1 measurement geometry for the spacer fluid are listed in Table 3. As in the previous example, the power law wall shear rate approximation produces shear rate estimates intermediate between the Newtonian and the Herschel-Bulkley estimates.  Table 1. We plot fits assuming both standard Newtonian shear rates, power law shear rates and Herschel-Bulkley shear rates.

Cement slurry
Finally, for the cement slurry, find very similar flow curve parametrizations, i.e. τ N = (1.07 + 1.10˙0 .69 ) Pa, for power law shear rates we obtain τ PL = (1.07 + 1.08˙0 .69 ) Pa, while utilizing the actual Herschel-Bulkley wall shear rates produces τ HB = (1.0 + 1.08˙0 .69 ) Pa. For completeness, we plot the Fann 35 measurements and the curve fits for the spacer fluid in Figure 9. The Newtonian, power law and Herschel-Bulkley shear rates in a Fann Model 35 R1-B1 measurement geometry for this cement slurry are listed in Table 4. Since this cement slurry is predominantly viscous with little degree of shear thinning and low yield stress value, the difference between Newtonian, power law and Herschel-Bulkley wall shear rates in the viscometer are small and the parametrizations virtually identical.  Figure 7 in the concentric annulus between a 9 5 ⁄ 8 " casing and a 12 1 ⁄ 4 " diameter wellbore. The pressure gradients are obtained for Herschel-Bulkley parametrizations based on different wall shear rate approximations.

Effect on friction pressure gradient; WBM
To gauge the significance of these fairly minor differences in flow curves, we evaluate the friction pressure gradient associated with fully developed laminar flow in the concentric annulus between a 9 5 ⁄ 8 " casing placed in a 12 1 ⁄ 4 " diameter wellbore, focusing on the water-based drilling fluid. These dimensions are common for production casings, one of the most important structural elements in wells today. Flow of this water-based drilling fluid in the annulus between casing and wellbore is relevant for e.g. preparation and execution of primary cementing operations; i.e. the operation where cement slurry is injected into the annulus displacing the original drilling fluid. We use the semi-analytical result of Hanks [18] to estimate the laminar friction pressure gradients and plot the results in Figure 10 for relevant flow rates. As expected, the difference in friction pressure gradient predictions is of the order of a few percent between the parametrizations obtained with Newtonian and Herschel-Bulkley wall shear rates; approximately 5.2% relative difference at a flow rate of 2000 l/min, with the Newtonian shear rate assumption producing the greatest friction pressure gradient. since the shear rate can be solved explicitly for a given combination of κ, Ω and n.

Discussion
Based on the examples shown in section 3, we have seen that increasing τy or decreasing k and n increase the difference between the frequently assumed Newtonian shear rate and the actual non-Newtonian shear rate. Increasing the yield stress and/or decreasing the shear thinning index results in flow curves that increasingly deviate from the linear Newtonian flow curve and increases the error made when assuming Newtonian behaviour. In Figure 11 we plot the wall shear rate ratio for the three fluids defined by the Fann 35 R1-B1 measurements in Table 1. Non-Newtonian effects on the wall shear rate are largest at lower rotational speeds and for the more shear thinning and higher yield stress fluids.
As illustrated by Figure 10, the Newtonian assumption results in an over-estimation of the actual friction pressure gradient. Since the actual friction pressure gradient is lower, this results in a lower circulation pressure in the annulus and a lower risk of hydraulic fracturing of the formation. In other operations, such as primary cementing under laminar flow conditions, the fluid friction pressure gradient is often used to design viscosity hierarchy among the fluids that are circulating and being displaced in the annulus outside the casing; the displacing fluid should be effectively more viscous (have higher friction pressure gradient) than the fluid it is displacing [4,5]. Due to e.g. narrow pressure margins toward formation fracturing pressure, the friction pressure gradients for successive fluids may be designed to be very close; depending on the non-Newtonian behaviour of the involved fluids, the result could be viscosity unstable situations leading to more intermixing and poorer displacement than planned for.
Finally, we remark that additional effects may influence the viscometric characterization of non-Newtonian fluids. Examples include thixotropic behavior and apparent wall slip, particularly at lower shear rates. Fluids with suspended weighting particles may in addition exhibit migration or sedimentation of solids during the measurement series which again may affect the measured viscosity of the fluid. These are examples of additional effects that can influence viscometer measurements of non-Newtonian fluids such as wellbore fluids.

Conclusion
In this paper we have considered the steady, laminar Couette (viscometric) flow of a Herschel-Bulkley fluid where the outer cylinder rotates. Following the derivation by Chatzimina et al. [15], we discuss the implicit relations between the rotational velocity of the outer cylinder, the cylinder diameter ratio, Herschel-Bulkley model parameters and the wall shear rate at the bob. We analyze Newtonian, power law and Herschel-Bulkley model wall shear rate predictions for three representative well construction fluids, and find that the common Newtonian shear rate assumption leads to under-estimation of the actual shear rate at the bob. The power law assumption results in wall shear rates intermediate between the Newtonian and the Herschel-Bulkley shear rates.
When used for viscosity characterization, the different wall shear rates result in different parameter combinations, as exemplified by the Fann 35 viscometer measurements presented in this paper. As the diameter ratio between the rotor and the bob in the industry-standard Fann Model 35 R1-B1 geometry is close to unity, deviations between the Newtonian, power law and Herschel-Bulkley wall shear rates are minor for the well construction fluids assessed here. For critical well operations, or when using a viscometer with a wider gap, correcting for non-Newtonian wall shear rates should be considered to obtain accurate viscosity characterization.