EXACT DETERMINATION OF THE THREE PARAMETERS

OF A HERSCHEL-BULKLEY FLUID

Wellington Campos
Geraldo Spinelli Ribeiro
PETROBRAS, CENPES/DIPLOT/SEPROT
Edimir Martins Brandão
PETROBRAS, CENPES/DIPLOT/SETEP
Ilha do Fundão, Q. 7, Rio de Janeiro, RJ, 21410-0900, Brazil

Abstract. This work presents an algorithm for the determination of the three parameters of a
Herschel-Bulkley fluid, namely, τo, n and K, for a given set of Fann readings. This problem is not
trivial, since the shear rate applied to the fluid in between the bob and rotor sleeve walls is a
function of the behavior index. The whole problem can be decomposed into two other problems,
which must be solved in alternation. The first problem consists of the determination of the shear
rate at the bob wall, for each value of the rotational speed, as a function of the three parameters
in the Herschel-Bulkley model. This is done through the use of Newton's method. The second
problem consists of the determination of the three Herschel-Bulkley parameters as a function of
the shear stresses and calculated shear rates at the bob wall by the use of both regression and
Newton's method. To start the iterative process, the fluid is assumed initially as Newtonian. The
whole problem converges quickly. The related problem of calculating the expected Fann
readings as a function of the Herschel-Bulkley parameters is also considered. Numerical
applications are presented and the error associated to the usual way of adopting Newtonian
shear rates is evaluated.

Key-words: Rheology, Herschel, Bulkley, Viscometer, Regression

1. INTRODUCTION

This study concerns the calculation of the three parameters of a Herschel-Bulkley fluid from
readings of Fann or Couette type viscometers. These viscometers comprise a stationary inner
cylinder and a rotational outer sleeve (Fig. 1). The inner cylinder is called the bob. As the outer
sleeve rotates, torque is transmitted to the bob by the fluid contained in the annular space between
the bob and the sleeve. The bob actuates a spring, bringing about a deflection of a needle by some
angle, which can be read on the dial reading. Figures 2 and 3 presents a more schematic
representation of this viscometer, and are going to be used in the developments of the theory.
 Figure 1 – Couette type viscometer. The bob is kept stationary while
the sleeve rotates. Torque is transmitted to the bob wall
by the fluid contained in the annulus.

Although this study could as well be extended to Searle type viscometers, in which the bob
rotates, rather than the outer sleeve, this will not be performed here, and is left for future works.
The fluid placed in the annular space, in between the bob and the sleeve, is assumed to
follow the Herschel-Bulkley model, also known as the yield-power-law model

τ = τ 0 + Kγ n (1)

where, τ is shear stress, γ is the shear rate and τ0, K and n are the three parameters of the
Herschel-Bulkley model which are called the yield point, the consistency index and the behavior
index, respectively.
Other rheological models can are used to relate τ to γ , e.g., the Bingham

τ = τ 0 + µ p γ (2)

where τ0 is the yield point and µp is the plastic viscosity; the Ostwald-de Waele or power-law

τ = Kγ n (3)

where K is the consistency index and n is the behavior index.

It should be emphasized that, although the rheological parameters share different names in
different models, their numerical values are not the same. This way, the yield point in the
Herschel-Bulkley model and the yield point in the Bingham model have different values for the
same fluid.
The Herschel-Bulkley model, the Bingham and the Ostwald-de Waele models, among others,
are fitted to the viscometer data through statistical regression. Hemphill et al. (1993) considered
the Herschel-Bulkley model more accurate than both the Bingham and the Power law for several
different fluid systems. That conclusion looks reasonable, since both the Bingham and the
Ostwald-de Waele models can be obtained from the Herschel-Bulkley model by setting n=1 or
τ0=0, respectively, in Eq. (1). The Bingham and the power law are two parameter models
whereas the Herschel-Bulkley is a three parameter model. There are other three parameter
models, as well, a very promising one being that presented by Robertson and Stiff (1976)

τ = A(γ + γ0 )n (4)

where A is a parameter related to consistency index, and γ 0 is a corrective parameter added to the
shear rate. The advantage of the Robertson-Stiff model lies in that it is more mathematically
tractable than the Herschel-Bulkley. The disadvantage in that the parameters have a more
cumbersome physical interpretation than the parameters in the Herschel-Bulkley model..
One problem posed for the Herschel-Bulkley model is how to calculate its three parameters,
the yield point, the consistency index and the behavior index, from readings of a rotational
viscometer. Other problem is how to calculate pressure drop through pipe and annular flows, in
both the laminar and the turbulent regimes, with these fluids. Another problem is how to
calculate heat transfer with these fluids. In this work only the first of these problems is
considered.
The literature contains a number of studies related with all the three problems referred above.
Hakki et al. (1992) have derived and expression for the axial laminar flow of a Robertson-Stiff
fluid in concentric annuli while stating that a similar expression could not be developed for a
Herschel-Bulkley fluid. Hemphill et al. (1993) used the Herschel-Bulkley model to study some
new drilling fluids, called mixed metal hydroxide, and found this model to be superior to both the
Bingham and Ostwald-de Waele fluids. Vinod et al. (1994) have conducted a very interesting
study of the flow of a Herschel-Bulkley fluid through an eccentric annuli, with co-existing
regions of laminar and turbulent regions. Bing et al. (1995) have used the Herschel-Bulkley
model to calculate surge and swab pressures during pipe trips in inclined and horizontal oil wells.
Merlo et al. (1995) have developed analytical expressions and a computer program for
calculating pressure drop in pipes and annulus using the Herschel-Bulkley model. Maglione et al.
(1996) attempted to use the standing pipe pressures to determine the Herschel-Bulkley
parameters at the rig site while drilling and oil well. The problem of finding the correct values of
the parameters of a Herschel-Bulkley fluid, nevertheless, have not been properly addressed in the
literature. The proper method for doing this is described in the sections ahead.

2. BACKGROUND

In this section , the fundamental equations and concepts related to the problem of finding
the three parameters of a Herschel-Bulkley fluid are presented. An implicit equation for the shear
rate and expressions for statistically fitting the Herschel-Bulkley model are developed. Both these
developments will need numerical techniques in their applications.
 Figure 2 – Lateral and top views of a Couette-type viscometer.

2.1 Shear stress and shear rate equation

The shear stress can be determined from the viscometer reading if the geometry of the
system and the value of the spring constant are known. The torque at a point in the annular space
between the bob and the sleeve is given by

T = τ rϕ (2πrL )r (5)

where τrφ is the shear stress, L is the bob length and r is the distance to the bob axis. The torque is
also related to the angle of deflection of the needle, θ , by

T = cθ (6)

where c is the spring constant. The spring constant c is chosen in a way as to make the viscometer
reading θ numerically equivalent to the viscosity of a Newtonian fluid, expressed in centipoise,
when the rotor sleeve is turned at a rotational speed of 300 rpm.
The combination of Eqs. (5) and (6) gives

cθ
τ rϕ = (7)
2πr 2 L

which shows that the shear stress varies inversely with the radius. It is possible to write
 τ r
2

τ rϕ = b 2b (8)
r

where rb is the bob radius and τb is the shear stress at the bob wall surface (see Fig. 2).
The velocity at a point in between the bob and the sleeve is given by

v = ωr (9)

where ω is the angular velocity and r is the distance to the bob axis. By taking the derivative

dv dω
= r +ω (10)
dr dr

it can be seem that the second term in the second hand side is related to translation, while the first
term in the second hand side is related to deformation, i.e., to the shear rate. Hence

dω
γ rϕ = r (11)
dr

Now, combining Eqs. (11) and (1) a differential equation is obtained

 dω 
n

τ rϕ = τ0 + Kr  (12)
 dr 

Integrating from the bob wall to the sleeve wall, using Eq. (8) and

2πN
ωs = (13)
60

where N is the rotational speed in rpm, the following expression ensues

1n
2πNK 1 n s 1  τ b rb 
r 2

= ∫  2 − τ 0  dr (14)
60 rb
r r 

An alternate expression is obtained by substituting Eq. (7) into (14)

2πNK 1 n s 1  cθ
r 1n

=∫  − τ 0  dr (15)
60 rb
r  2πr L
2


and even then an algorithm could be devised from this expression, the authors have taken a
different route. Equation ((15) is therefore, presented only for the sake of completeness.
It is convenient to change variables in Eq. (14)
 r
y= (16)
rb

so that it could be written

1n
2πNK 1 n 1 τb 
ys

60
= ∫1
 2 − τ 0  dy
yy 
(17)

Finally, substituting for the bob wall shear stress, the expression used in this work for the
calculation of the shear rate is derived

1n
1  τ 0 + Kγ b 
ys n

f (γ b ) =
60
2πNK 1 n ∫ y  y 2
− τ 0
 dy − 1 = 0

(18)
1 

This equation must be solved numerically, since the integral has no algebraic solution. In
addition, the expression inside the parenthesis cannot become negative, and therefore the
following condition must be fulfilled:

τ 0 ( y 2 − 1)
Π= ≤1 (19)
Kγ b
n

2.1 Regression equations

Expressions to fit the Herschel-Bulkley model to the viscometer data are developed in this
section. The viscometer data is a set of ordered pairs (Ni,θi), i=1...n, summarizing the viscometer
reading for each value of the rotational speed. Normally, n=6, with N taking values of 2, 3, 6,
100, 200, 300 and 600 rpm. In order to fit the Herschel-Bulkley model to the viscometer data, the
set (Ni,θi), i=1...n must be transformed into the set ( γi ,τbi), i=1...n, with the help of Eqs. (7) and
(18). Evidently, to perform this transformation, guessed values of τ0, K and n must be chosen.
The fitting of the Herschel-Bulkley equation requires the solution of the following minimization
problem

(
Min S (τ 0 , K , n ) = ∑ τ 0 + Kγi − τ bi
n
)2
(20)
i

where the index i is a label for each data point.

Taking advantage of the linear dependency of S in relation to τ0 and K, these parameters can
be written explicitly as
 ∑τ ∑ γ − ∑τ bi γ bi ∑ γ bi
2n n
bi bi
τ0 = i i i i
2
(21)
 n 
n ∑ γ bi −  ∑ γ bi 
2n

i  i 

n ∑τ bi γ bi − ∑ γ bi ∑τ
n n
bi

K=
i
2
(22)
 n 
n ∑ γ bi −  ∑ γ bi 
2n

i  i 

As for the behavior index, an implicit equation is obtained

g (n ) = τ 0 ∑ γi ln γ + K ∑ γi ln γ i − ∑τ bi γi ln γi = 0

2n n
(23)
i i i

4. NUMERICAL PROCEDURE

Given a data set consisting of viscometer readings θi at different rotational speed, Ni, i.e., (Ni,θi),
i = 1 to n, an algorithm can be formulated as a series of steps:

(a) Calculate the bob wall shear stresses, τbi, for each θi, using Eq. (7);
(b) Assume arbitrary values for τ0, K and n. Suggestion: n=1; τ0=0 N/m2; K=1 N sn/m2;
(c) Solve Eq. (18) using Gauss quadrature and Newton's method: get γi for each Ni;
(d) Solve equations (21) to (23) using Newton's method: get τ0, K and n;
(e) Compare this new values of τ0, K and n with the ones guessed in step (b);
(f) If the difference is too large, substitute the guessed values at step (b) by the calculated
values of step (d) and iterate from step (c);
(g) If the difference is sufficiently small, accept the values of τ0, K and n as the solution.

It may be convenient also to check the residue of Eqs. (18) and (23) after achieving
convergence in step (g). This algorithm should converge quickly if the data is legitimate. Of
course, it is always dangerous to input arbitrary data to any algorithm.
The authors have implemented this algorithm in two programming languages, FORTRAN
and PASCAL. They have also implemented the algorithm in Microsoft EXCEL, with the help of
the SOLVER capability.

5. RESULTS

The theory developed in the previous sections are now used to study a practical case. The
viscometer parameters are given in Table 1. The case examined consists in the calculation of the
Herschel-Bulkley parameters for an oil base mud studied by Hemphill et al. (1993). Table 2 and
Fig. 3 summarize the results. The residues of the functions f( γ i ), Eq. (18), and g(n), Eq. (23) are
shown close to zero as demanded by the theory. In the process of finding the zeros of these
functions, the condition for Π, expressed by Eq. (19), has been violated at the 3 rpm rotational
speed. A possible explanation is for the existence of a plug flow, or a region of no-shear in the
annular space, located close to the sleeve wall, where the shear stress is lower. This condition
would invalidate the theory developed previously.

Table 1 – Viscometer data.

Parameter Value
Spring constant, c 3,62615x10-5 N m/degree
Bob radius, rb 1,7245x10-2 m
Rotor sleeve radius, rs 1,8415x10-2 m
Bob length, L 3,8x10-2 m

Table 2 – Rheological parameters for the oil-base mud

through the present method of analysis.

Rotational Dial Bob Wall Bob Wall Function Parameter Herschel

Speed Reading Shear Rate Shear Stress Eq. (18) Eq. (19) Bulkley
Ni, rpm qi, degree γ i , sec-1 τbi, Pa f( γ i ) Πi Parameters
3 7 5,107 3,57 - 1,7429 n = 0.883
6 8 15,593 4,09 4,1x10-07 0,6504 K = 0.0628 Pa sn
100 18 178,942 9,19 1,5x10-05 0,0754 τ0 = 3.292 Pa
200 28 351,230 14,30 6,2x10-06 0,0416 g(n) = -5.2x10-7
300 38 523,172 19,41 -2,3x10-04 0,0292
600 63 1039,449 32,17 4,1x10-04 0,0159

35

30

25
Shear stress, Pa

20

15

10

0
0 200 400 600 800 1000 1200
Shear rate, 1/s

Figure 3 - Rheogram for the oil-based mud analyzed.

 Table 3 – Rheological parameters for the oil-base mud
through the classical method of analysis.

Rotational Dial Bob Wall Bob Wall Function Parameter Herschel

Speed Reading Shear Rate Shear Stress Eq. (18) Eq. (19) Bulkley
Ni, rpm qi, degree γ i , sec-1 τbi, Pa f( γ i ) Πi Parameters
3 7 5,107 3,57 - 1,7091 n = 0.876
6 8 10,214 4,09 2,3E-01 0,9311 K = 0.0666 Pa sn
100 18 170,229 9,19 2,7E-03 0,0791 τ0 = 3.388 Pa
200 28 340,459 14,30 1,0E-03 0,0431 g(n) = -1.3x10-8
300 38 510,688 19,41 6,4E-04 0,0302
600 63 1021,377 32,17 3,2E-04 0,0165

Consequently, the data for 3 rpm has been disregarded and the calculation proceeded for the
remaining five readings. The values of n, K and τ0 resulted somewhat different from those
obtained by Hemphill et al. (1993). This happened because Hemphill et al. (1993) have given up
the task of zeroing the two functions f and g, and stopped the iterations as soon as they found the
problem of evaluating the non-integer power of a negative value in Eq. (18).
Table 3 summarizes the results of the classical way of performing rheological calculations,
when the values of shear rate for a Newtonian fluid is used to fit the data of a non-Newtonian
fluid. It is known that, for Newtonian fluids, the shear rate do not depend on the absolute
viscosity. As can be shown in Table 3, the differences are not overwhelming, even though this
difference becomes lower and lower as the value of n tends to unity. In the case examined, the
relative errors are about 0,8%, 6% and 3% for the values of n, K and τ0 , respectively. These
errors are perfectly acceptable for many engineering applications. Nevertheless, it would be
interesting to investigate these same relative errors when the n values get farther from unity. As
Figure 3 demonstrate, the oil-base mud herein analyzed behaves almost as a Newtonian fluid.

5. CONCLUSIONS

A theory for determining the three parameter of a Herschel-Bulkley fluid has been developed
and an algorithm described to apply the theory. The authors have implemented the algorithm in
two programming languages, FORTRAN and PASCAL, and also in a spreadsheet, Microsoft
EXCEL. In this last case, the SOLVER capability has been used.
It has been shown through a case of the literature how to calculate the three parameters with
the equations herein developed. It is suggested that a plug flow ensues in the annular space close
to the sleeve wall at small shear rates. This violates the condition of validity of the theory. The
way around the problem was to disregard the data for small values of rotational speed and use the
remaining data to perform the calculations.
The comparison with the classical method of using Newtonian shear rates gave acceptable
errors. It remains to study the case when the value of n gets too far away from unity.

6. REFERENCES

Bing, Z., Kaiji, Z. and Qiji, Y., 1995, “Equations help calculate surge and swab pressures in
inclined wells,” Oil and Gas Journal, Sept. 1995, pp. 74-77.
Bingham, E. C., 1922, Fluidity and Plasticity, McGraw-Hill Book Co., New York, 1922.

Gücüyener, H. I. and Mehmetoglu, T., 1992, “Flow of Yield-Pseudo-Plastic Fluids through a

Concentric Annulus,” AIChE Journal, Vol. 38, No. 7, pp. 1139-1143.

Hemphill T., Campos, W. and Pilehvari, A., 1993, “Yield-power law model more accurately
predicts mud rheology,” Oil and Gas Journal, Aug. 23, pp. 45-50.

Maglione, R., Robotti, G. and Romagnoli, R., 1996, “A Computer Program to Predict Stand Pipe
Pressure while Drilling Using the Drilling well as Viscometer,” SPE 35994, presented at the
Eleventh Petroleum Computer Conference held in Dallas, Texas, June 2-5, 1996.

Robertson, R. E. and Still Jr., H. A., 1975, “An Improved Mathematical Model for relating Shear
Stress to Shear rate in Drilling Fluids and Cement Slurries,” SPE 5333, presented at the SPE-
AIME Rocky Mountain Regional Meeting held in Denver, April 7-9,1975.

Vinod, P. S., McIntire, L. V., “Rheology Effects on Laminar and Turbulent Flow of Model
Drilling Fluids Through Deviated Wells,” SPE 28287, provided to the Society of Petroleum
Engineers for distribution and possible publication in a SPE Journal, pp. 1-44.