energies

Article
Selection of a Suitable Rheological Model for Drilling
Fluid Using Applied Numerical Methods
Rafał Wiśniowski * , Krzysztof Skrzypaszek * and Tomasz Małachowski
Department of Drilling and Geoengineering, Faculty of Drilling Oil and Gas,
AGH University of Science and Technology, Mickiewicza 30 Av., 30–059 Kraków, Poland;
tomaszpmalachowski@gmail.com
* Correspondence: wisniows@agh.edu.pl (R.W.); varna@agh.edu.pl (K.S.)




Received: 8 May 2020; Accepted: 16 June 2020; Published: 19 June 2020


Abstract: The accuracy of fitting the rheological model to the properties of actual drilling fluid
minimises the errors of the calculated technological parameters applied while drilling oil wells.
This article presents the methodology of selecting the optimum drilling fluid rheological model.
Apart from classical rheological models, i.e., the Newtonian, Bingham Plastic, Casson, Ostwald
de Waele and Herschel–Bulkley models, it has been proposed to consider the Vom Berg and
Hahn-Eyring models, which have not been applied to describe drilling fluids so far. In the process
of determining rheological parameters for the Bingham Plastic, Casson, Ostwald de Waele and
Newtonian models, it is proposed to use a linear regression method. In the case of the Herschel–Bulkley,
Vom Berg and Hahn-Eyring models, it is suggested to use a non-linear regression method. Based on
theoretical considerations and mathematical relations developed in the Department of Drilling
and Geoengineering, Drilling, Oil and Gas Faculty, at AGH University of Science and Technology,
an original computer program called Rheosolution was developed, which enables automation of
the process of determining the optimum drilling fluid rheological model. Some examples show the
practical application of the method of selecting the optimum drilling fluid rheological model. Taking
into account data from actual measurements of drilling fluid properties, it has been proven that the
Vom Berg and Hahn-Eyring rheological models are best fitted to the description of drilling fluid
rheological parameters.

Keywords: drilling; drilling fluids; rheology; rheological model; numerical methods

1. Introduction
For over 150 years, different drilling fluids such as drilling muds, drill-in fluids, spacers and
washers, cement slurries and fracturing fluids have been used in the drilling practice. Each of them
has to fulfil strictly specified functions so that the borehole making process is safe both for the natural
environment and people working on a drilling rig. For instance, drilling mud is aimed at cleaning the
borehole bottom off drill cuttings, removing cuttings up to the area surface, ensuring the stability of
borehole walls both during drilling and circulation breaks, protecting against eruption (by exerting
hydrostatic pressure), maintaining cuttings in a state of suspension during breaks in circulation,
cleaning teeth of a drilling tool, lubricating drill bearings, cooling down a drill, and at supplying
energy to the borehole bottom in order to support the drilling process. In addition, during directional
drilling, horizontal drilling and multilateral drilling, drilling mud provides drive to mud motors
(PDM motors, turbine drilling motors, RSS systems) and enables control of the trajectory course
(Measurement-While-Drilling systems).
Together with the development of new drilling technologies such as high pressure and high
temperature (HPHT technology), underbalance drilling (UBD drilling), managed pressure drilling

Energies 2020, 13, 3192; doi:10.3390/en13123192 www.mdpi.com/journal/energies

Energies 2020, 13, 3192 2 of 17

(MPD technology), it becomes more and more important to determine precisely drilling mud and
drilling fluid parameters and properties. Technological properties of the applied fluids depend,
inter alia, on their rheological parameters. In order to describe cause and effect relations occurring
between rheological parameters of a fluid and the technology of its use, rheological models are
developed. The knowledge of a fluid rheological model is crucial since the accuracy of fitting a
rheological model to the actual property of fluid minimises errors of values being calculated, including
the nature of the fluid flow, fluid flow resistance in the circulatory system, diameters of drilling or
extending tool nozzles, mechanical and hydraulic parameters of the drilling technology, particle
sedimentation and the effectiveness of cuttings removal, injection radius and the technology of its
application (stream of flow volume, pressure and time of injection performance), power consumption
during mixing, mixing efficiency and issues relating to thermal conductivity during flow and mixing.
A number of rheological models have been developed in the drilling practice, however, none of
them is universal enough to describe precisely the behaviour of drilling fluids in a wide range of shear
rates. Rheological models form a mathematical description of the behaviour of Newtonian fluids.
Those models may be divided into two groups:
− relations between the shear rate and shear stress,
− relations between the shear rate or shear stress and apparent viscosity.
Models in both groups are different with respect to the applied function and parameters
approximating actual fluid behaviour. Models from the first group are applied for the analysis
of drilling fluids.

2. Rheological Models of Drilling Fluids

The following rheological models are most frequently applied for drilling technologies [1,2]:
!
dv
Newtonian : τ = η − (1)
dr
!
dv
Bingham : τ = τ y + ηpl − (2)
dr
!n
dv
Ostwald de Waele : (Power Law) : τ = k − (3)
dr
!1
2
1 1 1 dv
Casson : τ = τ y + ηcas
2 2 2
− (4)
dr
!n
dv
Herschel-Bulkley : (Yield-Power Law) : τ = τ y + k − . (5)
dr
The Bingham Plastic model has been applied most frequently due to its relative simplicity and
ease in calculating flow resistances [1]. Yet, at high shear rates and more complex drilling fluid
formulas, particularly those using polymers, this model does not describe the nature of drilling fluids
in a precise way.
When computers started to be used more commonly, it became possible to support computational
processes commonly and numerically, which helped to popularise the power model and the
Herschel—Bulkley model (YPL), which describes drilling fluids the best among those currently used
in the drilling practice [3,4]. Presently, computers are commonly available, and their computational
power enables the performance of quite complex and extended numerical calculations in a short time.
Therefore, it seems reasonable to analyse other models, which could prove to be better than those
applied so far, and the use of which could have been impractical so far, due to calculations complexity.
The Vom Berg model and the Hahn-Eyring model have been selected for the analysis [5,6].
The main selection criterion was another nature of the curve describing the dependency of shear stresses
Energies 2020, 13, 3192 3 of 17

on the shear rate. Contrary to the Ostwald—de Waele (Powel Law) model and the Herschel—Bulkley
(Yeld-Power Law) model, and also similar models, such as the Sisko model or the Robertson–Stiff
model [6], the curve in the analysed models is described by an inverse hyperbolic function instead
of a power function. The second criterion was very good results of applying those models in other
industries, which would justify considering them in the context of describing drilling fluids by means
of them. Both the Vom Berg model and the Hahn-Eyring model have successfully been applied to the
description of cement slurries used in the construction industry, which is confirmed in the literature [6].
The last condition was potentially high flexibility and possibility to fit into measurement data. In this
case (similarly to the YPL model), it is ensured by three parameters, which can be adapted to the fluid
under analysis.
New Vom Berg and Hahn-Eyring models, proposed to be applied in the drilling sector,
are described as follows [5–7]:
 dv 
 − 
Vom Berg : τ = τ y + Dsinh  dr 
−1 
 (6)
C
 dv 
 − 
!
dv
Hahn-Eyring : τ = E − + Dsinh  dr .
−1 
(7)
dr C
The authors make an attempt to demonstrate the usefulness of the Vom Berg and Hahn-Eyring
models (not used in calculations for drilling fluids yet) in minimising errors in the calculation of
rheological parameters and, in consequence, more accurate determination of flow resistance of
drilling fluids.

3. Methodology
Laboratory tests of real drilling fluids were performed by means of rotational viscometers.
The result of those tests was a set of measurement points, in which values of shear stresses (τ) occurring
. 
in the fluid under the influence of different shear rates (γ = − dr ) were determined.
dv

In order to determine the optimum rheological model, it has been suggested to conduct a
regression analysis, since it specifies the most probable function linking average values of both variables.
The procedure of determining a regression function is called the least-squares method. This principle
boils down to the statement that among all functions illustrating dispersion of measurement results the
best one is the one for which the sum of squares of measurement points deviations from that function
is the lowest [8].
In the case of the Bingham Plastic, Casson, Ostwald de Waele and Newtonian models, the authors
propose to apply a linear regression method. In case of the Herschel–Bulkley, Vom Berg and
Hahn-Eyring models, they suggest applying a non-linear regression method.

3.1. Linear Regression

In the case of linear regression, the optimum function is expected to take the form of a linear
Equation [8]:
ŷ = ax + b. (8)

The least squares estimation is presented in the following form:

m
X
U= ( yi − ŷ)2 → min . (9)
i=1
Energies 2020, 13, 3192 4 of 17

Taking into account Equation (8), a sum of squared residuals is obtained:

m 
X 
U= y2i − 2axi yi − 2byi + a2 x2i + 2abxi + b2 → min . (10)
i=1

Equation (10) is a function of coefficients a and b. The condition of function U minimisation can
be thus expressed in the form of the following equations:
 m m m
∂U
x2i + 2b
 P P P


 ∂a
= −2 xi yi + 2a xi = 0
i=1 i=1 i=1

. (11)

 m m
∂U
 P P



 ∂b
= −2 yi + 2a xi + 2bm = 0
i=1 i=1

By solving a system of Equation (11), one can acquire formulas for linear regression coefficients,
known from the literature: Pm Pm Pm
i=1 xi yi − i=1 xi i=1 yi
a= Pm P 2 (12)
2− m
i=1 x i i=1 x i
Pm Pm
yi − a i=1 xi
b = i=1 (13)
m
For the Bingham Plastic model, Equations (12) and (13) can be used directly, by applying
.
substitutions: η = a, τy = b, γ = x and τ= y. Rheological parameters are determined from the following
relation:
m
P m
P m
P
m xi yi − xi yi
i=1 i=1 i=1
η= !2 (14)
m m
2
P P
m xi − xi
i=1 i=1
m m
yi − η
P P
xi
i=1 i=1
τy = . (15)
m
For the Ostwald—de Waele (Power Law) model and the Casson model, their linearisation is
necessary. Linearisation of the Ostwald—de Waele (Power Law) model is acquired by logarithmic
Equation (3) on both sides as follows:
.
ln τ = ln k + n ln γ (16)
.
and using substitutions n = a, lnk = b, lnγ = x, lnτ = y.
Rheological parameters are determined from the following relation:
m
P m
P m
P
m xi yi − xi yi
i=1 i=1 i=1
n= !2 (17)
m m
x2i
P P
m − xi
i=1 i=1

m
P m
P
yi −η xi
( i=1 m i=1 )
k=e . (18)
Energies 2020, 13, 3192 5 of 17

√ √
Linearisations of the Casson model are performed by applying substitutions: ηcas = a, τ y = b,
q
. √
γ = x and τ = y. Rheological parameters are determined from the following relation:

 2
 P m m
P m
P 
 m xi yi − xi yi 
 i=1 i=1 i=1 
ηcas =  (19)
 
!2 
 m m 
 m x2i −
P P
xi
 
i=1 i=1

m m 2
√
P
η
P

 yi − cas x 
i

 i=1 i=1 
τ y =   . (20)


 m 


For the Newtonian model, when applying the least squares method, a regression function in the
following form is assumed:
ŷ = ax. (21)

The least-squares condition given with Equation (9), taking into account Equation (21) is described
with the formula:
Xm  
U= y2i − 2axi yi + a2 x2i → min . (22)
i=1

By differentiating Equation (22) against coefficient a, the following formula is obtained:

m m
∂U X X
= −2 xi yi + 2a x2i (23)
∂a
i=1 i=1

from which, after comparing to zero, one can designate the value of coefficient a
m
P
xi yi
i=1
a= m
. (24)
x2i
P
i=1

A rheological parameter of the Newtonian model (dynamic viscosity η) is thus determined from
the following relationship:
m . 
γi τi
P
i=1
η= m . (25)
P .2
γi
i=1

3.2. Non-Linear Regression

The Herschel–Bulkley model cannot be linearised since while determining parameter equations
an implicit equation with one variable occurs. The condition of function U minimisation can thus be
expressed in the form of equations:
 m m .n
∂U
τi + 2k γi + 2mτ y = 0
P P
= −2


∂τ y



 i=1 i=1
m m m

∂U .n .n . 2n

τi γi γi + 2k
γi = 0
 P P P

 ∂k
= −2 + 2τ y . (26)
i=1 i=1 i=1



 m .n . m
P .n . m
P . 2n .
∂U
    
= −2k (τi γi ln γi ) + 2τ y k γi ln γi + 2k2 γi ln γi = 0

 P
∂n


i=1 i=1 i=1
Energies 2020, 13, 3192 6 of 17

A system of Equation (26) after ordering variables boils down to the following form: [4]
m .n ..n m
. m
(τi γi ln γi ) −
(γi ln γi )+τi
P P P
m
i=1 i=1 i=1
m .n Pm m .n
m τi γi − τi
P P
γi " m
P .nP m .n m  . 2n
# (27)
.  . 
i=1 i=1 i=1
γ γ γ γ γ
P
m . 2n Pm .n Pm .n i i ln i − m i ln i = 0.
m γi − γi γi i=1
P
i=1 i=1
i=1 i=1 i=1

Equation (27) is an implicit equation with one variable n. In order to solve it, it is suggested to
apply one of the numerical methods, e.g., the bisection method [9].
In this method, a function is constructed in the following form: [4]
m .n . m m .n . 
(τi γi ln γi ) − τi γi ln γi +
P P P
g(n) = m
i=1 i=1 i=1
m .nm .n m
m
P
τi γi −
P
γi " m
P
τi m
P .nP .n .  m  . 2n
# (28)
 . 
i=1 i=1 i=1
γ γ γ γ γ
P
m . 2n P
m .n P
m .n i i ln i − m i ln i
m γi − γi γi i = 1
P
i=1 i=1
i=1 i=1 i=1

and next zero is designated.

To this end, at the beginning, the range [A, B] is established, within which the existence of the root
being searched for is predicted. Considering the physical sense of the parameter n being searched for,
one can assume n ∈ [0, 100]. The principle of designating the zero of the g(n) function in the range
[A, B] by means of the bisection method is as follows (Figure 1):
Energies 2020, 13, x FOR PEER REVIEW 7 of 20

Figure 1. Determination
Figure 1. Determination of
of the
the root
root of the g(n)
of the g(n) function
function by
by means
means of
of the
the bisection
bisectionmethod.
method.

The starting range [A, B] is divided into two in order to determine the value of the coordinate
The starting range [A, B]+isB divided into two in order to determine the value of the coordinate of
of the interval centre: nsr =AA+ 2 B. Next, the value of the g(nsr ) function is determined and compared
the interval centre: nsr = . Next, the value of the g(nsr) function is determined and compared
to the g(A) value. If the function 2 changes the sign, it means that the root of the g(n) function is in
the range [A, n ]; otherwise, the root
śr If the function changes
to the g(A) value. is inthe
thesign,
range [nśr , B].
it means Calculations
that the root of are
the made until theis range
g(n) function in the
width
range [A, nśr]; otherwise, the root is in the range [nśr, B]. Calculations are made until the rangeiswidth
is smaller or equal to accuracy ε assumed at the beginning. The root of the function then
is smaller or equal to accuracy ε assumed at the beginning. The root of the function is then determined
based on the formula: nsr = A + B , whereas the error of its estimation is calculated from the
2
dependency: blad = B − A .
2
Energies 2020, 13, 3192 7 of 17

determined based on the formula: nsr = A+ B

2 , whereas the error of its estimation is calculated from the
B−A
dependency: blad = 2 .
The algorithm of Equation (28), with the use of the bisection method, has been provided in
Energies 2020,
Figure 2. 13, x FOR PEER REVIEW 8 of 20

g(x)

Start
equation 28
Input data
[i],[i],m,A,B, epsilon
return

nsr=(A+B)/2

gsr=g(nsr)
ga=g(A)

No No
ga*gsr=0 ga*gsr<0 A=nsr
Yes Yes
No
ga=0 B=nsr
Yes

nsr=A No
nsr=b
|B-A|=epsilon

Yes

nsr=(A+B)/2

ERROR=0 ERROR=(B-A)/2

Print
nsr, ERROR

Stop

Figure
Figure 2. Flow chart
2. Flow chart for
for the
the calculation
calculation of
of the
the root
root of the g(n)
of the g(n) function
function within the interval
within the [A, B]
interval [A, B] by
by
means of the bisection method.
means of the bisection method.

Given the designated value of the parameter n, the other rheological parameters in the
Given the designated value of the parameter n, the other rheological parameters in the Herschel-
Herschel-Bulkley model can be specified from Equations (29) and (30):
Bulkley model can be specified from Equations (29) and (30):
m
m .n mm mm .n
 − i
m τiγin −
1 i
τi γi n
P P P
m i=1 i i =1 ii=
kk == i =m1
P . 2n P mi =1. n P
mi =1 . n
(29)
(29)
m m γi2 n − m γi n m γi n
m
i=1 i − i 1 i
=1 i i=
i =1 i =1 i =1
m m.n
τi − k γi
P P
m m

 − k  
i=1 i=1 n
τy = . (30)
i m i
(30)
 y =model, also in the. case of the Vom Berg and Hahn-Eyring
Similarly, to the Herschel-Bulkley (YPL)
i =1 i =1
m
models, linearisation of their rheological equations is not possible:
The condition
Similarly, of the
to the function U minimisation
Herschel-Bulkley for the
(YPL) model, alsoVom Berg
in the model,
case using
of the Vom substitutions
Berg and Hahn-
Eyring models, linearisation of their rheological equations is not possible: !
The condition of the function U minimisation for the Vom Bergdvmodel, using substitutions
a = τ y , b = D, c = C, yi = τi , xi = −
𝑑𝑣 dr i
𝑎 = 𝜏𝑦 , 𝑏 = 𝐷, 𝑐 = 𝐶, 𝑦𝑖 = 𝜏𝑖 , 𝑥𝑖 = −( )𝑖
𝑑𝑟
can be expressed in the form of a system of equations:
Energies 2020, 13, 3192 8 of 17

can be expressed in the form of a system of equations:

∂U P   x 
=2 m bsinh−1 ci − 2 m
 P
∂a i = 1 i=1 yi +2am = 0



∂U Pm  −1 xi −1 xi
   
= i=1 −2sinh c −a − bsinh c + yi = 0

∂b


. (31)

−1 xi

 ∂U P m 2bx i ( −a−bsinh ( c ) + y i )


 ∂c
= i=1
r
2
= 0

 x
 c2 1+ i
c2

The condition of the function U minimisation for the Hahn-Eyring model, by applying substitutions
!
dv
a = E, b = D, c = C, yi = τi , xi = −
dr i

can be expressed in the form of a system of equations:

∂U −1 xi
   
= m
 P
∂a
( −2x i −ax i − bsinh + yi ) =0
Pm  i=1 c  



∂U −1 xi −1 xi
   
= i=1 −2sinh c −axi − bsinh c + yi = 0

∂b


. (32)


 ∂U Pm 2bxi (−axi −bsinh−1 ( xci )+ yi )


 ∂c
= i=1
r
2
=0

 x
 c2 1+ i
c2

There are a number of available non-linear regression methods [10,11]. With regard to the Vom
Berg and Hahn-Eyring models, it is proposed to apply gradient descent [12]. This algorithm enables us
to find local minima of a function and to determine the global minimum from among them. It consists
of selecting a starting point to be found in the domain of the analysed function. This function is a
function which specifies fitting of the analysed model equation (the sum of squared residuals) and
designates a vector in it with the direction of a given function gradient and the sense opposite to it.
The length of a step (vector) is equal to one multiplied by the vector length coefficient. This coefficient
can be freely adapted to optimise the algorithm:
n
X
U (a, b, c) = ( f (a, b, c) − yi )2 → min (33)
i=1

→
− ∇U (a, b, c)
v̂ = →
·αk . (34)
− ∇U (a, b, c)

After the computational sequence has been performed, the model parameters are modified with
the vector value and the algorithm selects a given point as a new starting point. The process is repeated
until a local minimum is achieved, with the assumed precision. Then, the starting point changes and
another local minimum is searched for. After repeating calculations, the assumed number of times,
the smallest of the local minima being found is selected and considered the global minimum.
The flow chart of calculations is presented in Figures 3 and 4.
To simplify the equation, particular rheological model parameters are ascribed to variables a, b, c,
whereas constant yi and xi are data from subsequent measurement points. With every step, only the
variables are modified. Assignment of model equations to the diagram:
 dv 
 − 
the Vom Berg model : τ = τ y + Dsinh  dr  −1 

C
!
dv
a = τ y , b = D, c = C, yi = τi , xi = −
dr i
Energies 2020, 13, 3192 9 of 17

m m  2
xi
X X  
U (a, b, c) = ( f (a, b, c) − yi )2 = a + bsinh−1 − yi → min (35)
c
i=1 i=1
 dv 
 − 
!
dv
the Hahn-Eyring model : τ = E − + Dsinh  dr 
−1 
dr C
!
dv
a = E, b = D, c = C, yi = τi , xi = −
dr i
m m  2
xi
X X  
U (a, b, c) = ( f (a, b, c) − yi )2 = axi + bsinh−1 − yi → min (36)
c
i=1 i=1

Partial derivatives are calculated according to Equation (31) or Equation (32). The gradient descent
can be also applied
Energies 2020, 13, xwhen designating
FOR PEER REVIEW rheological parameters under the Herschel-Bulkley model.
10 of 20

Figure
Figure 3. Flow
3. Flow chartdetermining
chart determining the
thefunction minimum.
function minimum.
Energies 2020, 13, 3192 10 of 17
Energies 2020, 13, x FOR PEER REVIEW 11 of 20

Figure 4. Schematic diagram of the non-linear programming method.

Figure 4. Schematic diagram of the non-linear programming method.

3.3. Numerical Support of the Process of Selecting a Drilling Fluid Rheological Model
To simplify the equation, particular rheological model parameters are ascribed to variables
𝑎, 𝑏, 𝑐, whereasofconstant
Development 𝑦𝑖 and 𝑥𝑖 are
the methodology ofdata fromrheological
fitting subsequent measurement
parameters in points. With every
research data step,
has enabled
only the variables are modified. Assignment of model equations to the diagram:
us to create a numerical tool making it possible to select the optimum drilling fluid rheological model.
𝑑𝑣
A new version of the Rheosolution program,
the Vom version
Berg model: 𝜏 =4.0,
𝜏𝑦 +has been
𝐷 sinh −1 developed,
−
( 𝑑𝑟) in which apart from the
𝐶
Newtonian, Bingham Plastic, Ostwald–de Waele, Casson and Herschel–Bulkley models, handling of
the Vom Berg model and the Hahn-Eyring model was introduced. 𝑑𝑣 A number of other changes were
𝑎 = 𝜏𝑦 , 𝑏 = 𝐷, 𝑐 = 𝐶, 𝑦𝑖 = 𝜏𝑖 , 𝑥𝑖 = −( )𝑖
made, including, inter alia, modernisation and reconstruction, 𝑑𝑟 together with rewriting the program
𝑚 𝑚
from the Pascal language into C++/CLI and also adding the English 𝑥𝑖 version.
2
𝑈(𝑎, 𝑏, 𝑐) = ∑(𝑓(𝑎, 𝑏, 𝑐) − 𝑦𝑖 = ∑ (𝑎 + 𝑏 sinh−1 ( ) − 𝑦𝑖 ) → 𝑚𝑖𝑛
)2
(35)
Rheosolution in the new version 𝑖=1
4.0 adapts the 𝑖=1
parameters 𝑐of seven rheological models to the
measurement data being introduced. The program calculates rheological parameters for each of the
Energies 2020, 13, x FOR PEER REVIEW 𝑑𝑣 12 of 20
𝑑𝑣 degree−1 −
analysed models and indicates the modelmodel:
the Hahn-Eyring with the𝜏 = highest
𝐸 (− ) + 𝐷 sinh of( fitting.
𝑑𝑟
) Data and diagrams are
𝑑𝑟 𝐶
exported to files,models
analysed images andand spreadsheets.
indicates the model with the highest degree of fitting. Data and diagrams are
The way oftofitting
exported rheological
files, images parameters depends on a model.
and spreadsheets. 𝑑𝑣 In case of the Newtonian and
The way of fitting 𝑎 = 𝐸,parameters
rheological 𝑏 = 𝐷, 𝑐 = 𝐶, 𝑦𝑖 = 𝜏𝑖 ,on
depends 𝑥𝑖 = −( )𝑖 In case of the Newtonian and
a model.
Bingham Plastic models, linear regression is applied, which enables
𝑑𝑟 designation of linear equation
Bingham Plastic models, linear regression is applied, which enables designation of linear equation
parameters. Ostwald–de Waele 𝑚 and Casson equations 𝑚 are linearised, and 2
next, linear regression is
parameters. Ostwald–de Waele and Casson 2equations are linearised, −1
𝑥𝑖 and next, linear regression is
applied. The last three 𝑈(𝑎,models
𝑏, 𝑐) = ∑(𝑓(𝑎, 𝑏, 𝑐) − 𝑦𝑖 ) by
are designated =∑ means
(𝑎𝑥 + of𝑏 gradient
sinh ( )descent.
− 𝑦𝑖 ) → 𝑚𝑖𝑛 (36)
applied. The last three models are designated by means𝑖 of gradient descent.
𝑐
The user interface has been𝑖=1 provided in Figure 𝑖=1 5. The method of program operation is shown in
The user interface has been provided in Figure 5. The method of program operation is shown in
Figure 6. Partial
Figure 6. derivatives are calculated according to Equation (31) or Equation (32). The gradient
descent can be also applied when designating rheological parameters under the Herschel-Bulkley
model.

3.3. Numerical Support of the Process of Selecting a Drilling Fluid Rheological Model
Development of the methodology of fitting rheological parameters in research data has enabled
us to create a numerical tool making it possible to select the optimum drilling fluid rheological model.
A new version of the Rheosolution program, version 4.0, has been developed, in which apart from
the Newtonian, Bingham Plastic, Ostwald–de Waele, Casson and Herschel–Bulkley models, handling
of the Vom Berg model and the Hahn-Eyring model was introduced. A number of other changes
were made, including, inter alia, modernisation and reconstruction, together with rewriting the
program from the Pascal language into C++/CLI and also adding the English version.
Rheosolution in the new version 4.0 adapts the parameters of seven rheological models to the
measurement data being introduced. The program calculates rheological parameters for each of the

Figure 5. Rheosolution 4.0 interface.

Figure 5. Rheosolution 4.0 interface.
Energies 2020, 13, 3192 11 of 17
Energies 2020, 13, x FOR PEER REVIEW 13 of 20

Figure 6.
Figure 6. Rheosolution
Rheosolution 4.0
4.0block
blockdiagram.
diagram.

4. 4.
Results
Resultsand
andDiscussion
Discussion
InIn
order
ordertotoexemplify
exemplifythetherelations
relations presented inthis
presented in thisarticle,
article,laboratory
laboratory tests
tests of of different
different drilling
drilling
fluids were made [1]. The measurements were made in the Drilling Fluids Laboratory
fluids were made [1]. The measurements were made in the Drilling Fluids Laboratory in the Department in the
of Department
Drilling andofGeoengineering,
Drilling and Geoengineering, Faculty of
Faculty of Drilling, OilDrilling,
and Gas Oilon
and Gas on
a Fann a Fann viscometer.
viscometer. The tested
The tested
fluids fluids
were the were the following:
following:
• • Cementslurry
Cement w/c==0.5,
slurryoow/c 0.5, without
without any additives,
additives,
• •Cement slurry o w/c = 0.5 with addition of
Cement slurry o w/c = 0.5 with addition of 0.3%
0.3%PSP
PSP042
042superplasticizer,
superplasticizer,
• • 3% bentonite drilling mud without any additives.
3% bentonite drilling mud without any additives.
• •3%3% bentonitedrilling
bentonite drillingmud
mudwith
withan
an addition
addition of
of 2%
2% XCD
XCDpolymer.
polymer.
Table 1 shows the results of rheological measurements of four drilling fluid formulas performed
Table 1 shows the results of rheological measurements of four drilling fluid formulas performed
on a Fann viscometer, with the arrangement of Rotor-Bob cylinders (R1-B1) and an applied spring
on a Fann viscometer, with the arrangement of Rotor-Bob cylinders (R1-B1) and an applied spring
Energies 2020, 13, 3192 12 of 17

(F1) [1]. Table 2 shows relations between shear stresses and the shear rate obtained from laboratory
tests on a Fann viscometer, with the arrangement of Rotor-Bob cylinders (R1-B1) and an applied spring
(F1) [1]. Tables 3 and 4 present the calculated values of rheological parameters pertaining to fluids
under analysis, approximated with particular rheological models. Tables 5 and 6 present rheological
model selection results, pertaining to the analysed fluids, obtained from Rheosolution. Figures 7–10
present drawings to show a comparison of the selected rheological models for the tested drilling mud
and cement slurry. To maintain the clarity of the drawing, models recommended by the American
Petroleum Institute (Bingham Plastic and Ostwald de Waele models) were selected for comparison.

Table 1. Data from laboratory tests made on a Fann viscometer (R1-B1, F1).

Cement w/c 0.5 with Addition of PSP 042 3% Drilling Mud with Addition of 2% XCD
Rotations Angle Rotations Angle Rotations Angle Rotations Angle
rot/min ◦ rot/min ◦ rot/min ◦ rot/min ◦

1 12 1 2 1 3 1 23
2 15 2 2 2 4 2 24
3 18 3 3 3 5 3 26
6 25 6 4 6 6 6 28
10 32 10 6 10 6 10 31
20 47 20 10 20 6 20 35
30 54 30 14 30 6 30 38
60 67 60 28 60 8 60 45
100 80 100 48 100 10 100 51
200 108 200 91 200 14 200 65
300 132 300 141 300 18 300 76
600 183 600 237 600 28 600 98

Table 2. Relations between shear stresses and shear rates obtained from laboratory tests on a Fann
viscometer (R1-B1, F1).

Cement w/c 0.5 with Addition of PSP 042 3% Drilling Mud with Addition of 2% XCD
Shear Rate Tension Shear Rate Tension Shear Rate Tension Shear Rate Tension
1/s Pa 1/s Pa 1/s Pa 1/s Pa
1.7034 6.132 1.7034 1.022 1.7034 1.533 1.7034 11.753
3.4068 7.665 3.4068 1.022 3.4068 2.044 3.4068 12.264
5.1102 9.198 5.1102 1.533 5.1102 2.555 5.1102 13.286
10.2204 12.775 10.2204 2.044 10.2204 3.066 10.2204 14.308
17.034 16.352 17.034 3.066 17.034 3.066 17.034 15.841
34.068 24.017 34.068 5.11 34.068 3.066 34.068 17.885
51.102 27.594 51.102 7.154 51.102 3.066 51.102 19.418
102.204 34.237 102.204 14.308 102.204 4.088 102.204 22.995
170.34 40.88 170.34 24.528 170.34 5.11 170.34 26.061
340.68 55.188 340.68 46.501 340.68 7.154 340.68 33.215
511.02 67.452 511.02 72.051 511.02 9.198 511.02 38.836
1022.04 93.513 1022.04 121.107 1022.04 14.308 1022.04 50.078

Table 3. Fluid rheological parameters obtained from Rheosolution 4.0 for cement slurry samples (with
and without adding superplasticizer).

Cement w/c 0.5 with Addition of PSP 042

Newtonian η = 0.1105 - - η = 0.1246 - -
Bingham ηpl = 0.0842 τy = 16.9723 - ηpl = 0.1215 τy = 1.9507 -
Ostwald—de Waele k = 4.7932 n = 0.4265 - k = 0.3931 n = 0.8012 -
Casson ηcas = 0.0532 τy = 8.7314 - ηcas = 0.1159 τy = 0.1668 -
Herschel- Bulkley k = 2.2527 n = 0.5420 τy = 1.2854 k = 0.2925 n = 0.8718 τy = 0.1440
Vom Berg τy = 9.8527 D = 26.5833 C = 105.2867 τy = 0.0165 D = 96.9329 C = 643.4336
Hahn-Eyring E = 0.0922 D = 1.6225 C = 0.7939 E = 0.1146 D = 0.9401 C = 1.0184
Energies 2020, 13, 3192 13 of 17

Table 4. Fluid rheological parameters obtained from Rheosolution 4.0 for drilling mud samples (with
and without adding polymer).

3% Drilling Mud with Addition of 2% XCD

Newtonian η = 0.0159 - - η = 0.0622 - -
Bingham ηpl = 0.0119 τy = 2.5903 - ηpl = 0.0377 τy = 15.8508 -
Ostwald—de Waele k = 1.3046 n = 0.2959 - k = 8.9200 n = 0.2241 -
Casson ηcas = 0.0057 τy = 1.7405 - ηcas = 0.0145 τy = 11.8160 -
Herschel-Bulkley k = 0.1682 n = 0.6217 τy = 1.3479 k = 4.6678 n = 0.3285 τy = 2.8143
Vom Berg τy = 0.4649 D = 1.8184 C = 9.2868 τy = 12.6203 D = 11.5721 C = 98.7531
Hahn-Eyring E = 0.0086 D = 0.6916 C = 1.0548 E = 0.0244 D = 2.4781 C = 0.0507

Table 5. Results of the rheological model selection obtained from Rheosolution 4.0 for drilling fluid
samples (with and without adding polymer).

3% Drilling Mud with Addition of 2% XCD

Model/Coefficient R F U R F U
Newtonian 0.77745 15.27968 59.73121 0 0 2256.306
Bingham 0.99074 532.5236 2.78326 0.96043 118.8995 123.9423
Ostwald—de Waele 0.92387 58.27622 22.11583 0.97087 164.1754 91.72426
Casson 0.99554 1113.56 1.34393 0.99417 849.8189 18.58076
Herschel-Bulkley 0.98992 488.5667 3.02866 0.98816 414.6451 37.62234
Vom Berg 0.88438 35.90033 32.89705 0.99375 792.1389 19.91689
Hahn-Eyring 0.98781 402.6485 3.65926 0.99679 1547.903 10.25482

Table 6. Results of the rheological model selection obtained from Rheosolution 4.0 for cement slurry
samples (with and without adding superplasticizer).

Cement w/c 0.5 with Addition of PSP 042

Model/Coefficient R F U R F U
Newtonian 0.77729 15.26429 3214.332 0.99534 1066.001 143.4055
Bingham 0.95144 95.52491 769.5607 0.99639 1378.773 111.1078
Ostwald—de Waele 0.99855 3449.434 23.47425 0.98011 243.846 607.8673
Casson 0.97705 210.3658 368.5144 0.99716 1750.87 87.62985
Herschel- Bulkley 0.98691 374.3877 211.2646 0.99874 3975.655 38.71535
Vom Berg 0.99009 496.8507 160.2201 0.99956 11232.52 13.7246
Hahn-Eyring 0.90418 44.80884 1481.655 0.99692 1616.03 94.89674

The tests were successful. In the case of pure slurry, the Ostwald–de Waele model turned out to
be the best-fitted model. For slurry with addition of PSP 042 superplasticizer, the Vom Berg model was
the best fitted model. Pure bentonite drilling mud was best described by means of the Casson model,
and polymer inhibited mud, by means of the Hahn-Eyring model. For drilling muds and cement
slurries without additives, differences between the correlation coefficient are not that significant among
some models and linear models have a good correlation in these cases.
Preliminary studies on sample drilling mud and cement slurries encourage more accurate
analysis; the results indicate that the multi-parameter models of Vom Berg and Hahn-Eyring have
a high correlation for drilling fluids with an addition of polymers or superplasticizers. This is due
to the fact that such additions most often cause non-linear increases in the value of rheological
parameters, therefore, models based on exponential functions better illustrate these changes. The use
of computational methodology, Rheosolution 4.0 and viscometers such as Brookfield R/S (multi shear
rate measurement points) can improve the adjustment of a rheological model to the behaviour of the
actual drilling fluid during its flow.
Energies 2020, 13,
Energies 2020, 13, 3192
x FOR PEER REVIEW 14 of
16 of 17
20
Energies 2020, 13, x FOR PEER REVIEW 16 of 20

Figure 7.
Figure 7. Comparison
Comparison of
of the
the selected
selected rheological models, cement
rheological models, cement w/c
w/c 0.5.
0.5.
Figure 7. Comparison of the selected rheological models, cement w/c 0.5.

Figure Comparison
8. 8.
Figure Comparisonofofthe
theselected
selectedrheological
rheological models, cement w/c
models, cement w/c0.5
0.5with
withan
anaddition
additionofofPSP
PSP 042.
042.
Figure 8. Comparison of the selected rheological models, cement w/c 0.5 with an addition of PSP 042.
Recently, progress has been made in the use of nanoparticles as an addition to drilling fluids in
order to adjust their rheological parameters. Authors [13,14] emphasise in their works that the Power
Law model or the Herschel Bulkley model can be used to predict drilling fluids behaviour during flow,
whereas the application of the Bingham Plastic model can sometimes only capture the trend. This is a
group of new drilling fluids which might require further studies to optimise the best mathematical
model. The issue seems to be interesting and worth analysing since these drilling fluids, as indicated
Energies 2020, 13, 3192 15 of 17

in some papers [15,16], describe exponential models well and such rheological models include also the
Energies 2020, 13, x FOR PEER REVIEW 17 of 20
Vom Berg
Energies and
2020, 13,Hahn-Eyring models.
x FOR PEER REVIEW 17 of 20

Figure9.9.Comparison
Figure Comparisonofofthe
theselected
selectedrheological
rheologicalmodels,
models,3%
3%drilling
drillingmud.
mud.
Figure 9. Comparison of the selected rheological models, 3% drilling mud.

Figure 10. Comparison of the selected rheological models, 3% drilling mud with an addition of 2% XCD.
Figure 10. Comparison of the selected rheological models, 3% drilling mud with an addition of 2%
Figure 10. Comparison of the selected rheological models, 3% drilling mud with an addition of 2%
XCD.
5. Conclusions
XCD.
For
Themany
testsyears, different drilling
were successful. In thefluids
case ofsuch asslurry,
pure drillingthe
muds, drill-in fluids,
Ostwald–de Waelespacers and washes,
model turned out to
cement The tests were
slurries and successful.fluids
fracturing In thehave
case been
of pure slurry,
applied inthe
theOstwald–de
drilling WaeleEach
practice. model turned
kind out to
of drilling
be the best-fitted model. For slurry with addition of PSP 042 superplasticizer, the Vom Berg model
be the
fluid best-fitted model. with
For slurry with addition of PSP 042 superplasticizer, ofthe Vom Berg model
was should
the bestbefitted
described
model. Pure thebentonite
best fitted rheological
drilling mud model.
was bestThe accuracyby
described fitting
means aofrheological
the Casson
was the best fitted model. Pure bentonite drilling mud was best described by means of the Casson
model, and polymer inhibited mud, by means of the Hahn-Eyring model. For drilling muds and
model, and polymer inhibited mud, by means of the Hahn-Eyring model. For drilling muds and
cement slurries without additives, differences between the correlation coefficient are not that
cement slurries without additives, differences between the correlation coefficient are not that
significant among some models and linear models have a good correlation in these cases.
significant among some models and linear models have a good correlation in these cases.
Energies 2020, 13, 3192 16 of 17

model to the properties of actual fluid minimises errors in values being calculated, such as the nature
of fluid flow, fluid flow resistances in the circulatory system, diameters of drilling or extending tool
nozzles, mechanical and hydraulic parameters of the drilling technology, particle sedimentation and
the effectiveness of cuttings removal, injection radius and the technology of its application (stream of
flow volume, pressure and time of injection performance), power consumption during mixing, mixing
efficiency and issues relating to thermal conductivity during flow and mixing.
The methodology developed and presented in this article enables us to select the optimum
rheological model of any drilling fluid. In the process of designating rheological parameters for the
Bingham Plastic, Casson, Ostwald de Waele and Newtonian models, it is proposed to apply a linear
regression method. In the case of the Herschel–Bulkley, Vom Berg and Hahn-Eyring models, it is
suggested to apply a non-linear regression method.
An interesting issue seems to be the use of the Vom Berg and the Hahn-Eyring rheological models
to describe the flow of drilling fluids with an addition of nanoparticles [13–16] and this will be one of
the later stages in the work of the authors of this article. We draw such conclusions on the basis of
observations that multi-parameter models better describe the flow of liquids modified with chemical
additives (superplasticizers, polymers).
The original Rheosolution computer program, developed in the Department of Drilling and
Geoengineering, Drilling, Oil and Gas Faculty, at AGH University of Science and Technology, enables
automation of the process of designating the optimum drilling fluid rheological model. Owing to
the mathematical analysis conducted in this study and numerical calculations, taking into account
data from actual measurements of drilling fluid properties, it has been proven that the Vom Berg
and Hahn-Eyring rheological models fit best to describe drilling fluid rheological parameters. It is
worthwhile to take subsequent steps in order to further test those models in terms of their usefulness
in the drilling industry.

Author Contributions: Conceptualization, R.W. and K.S.; methodology, R.W., K.S. and T.M.; software, T.M. and
K.S.; validation, R.W., K.S. and T.M.; formal analysis, R.W., K.S. and T.M.; investigation, R.W., K.S. and T.M.;
resources, R.W., K.S. and T.M.; data curation, R.W., K.S. and T.M.; writing—original draft preparation, R.W. and
K.S.; writing—review and editing, R.W. and K.S.; visualization, R.W., K.S. and T.M.; supervision, R.W. and K.S.;
project administration, R.W. and K.S. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.
Conflicts of Interest: The author declares no conflict of interest.

Nomenclature
a, b—regression model coefficients; [-]
c—regression model coefficient for the Vom Berg and Hahn-Eyring models; [-]
A, B—interval boundaries in the bisection method; [-]
C—rheological parameter in the Vom Berg and Hahn-Eyring models; [s−1 ]
D—rheological parameter in the Vom Berg and Hahn-Eyring models; [Pa]
E—rheological parameter in the Hahn-Eyring model; [P·s]
F—Fisher-Snedecor index; [-]
U—sum of squared residuals; [(Pa)2 ]
.
dv/dr—shear rate gradient (γ); [s−1 ]
ηpl —plastic viscosity; [Pa·s]
ηcas —viscosity of Casson; [Pa·s]
.
γi —shear rate measured at i-th rotational speed; [s−1 ]
k—coefficient of drilling mud consistency; [Pa·sn ]
m—number of measurements with viscometer; [-]
n—exponential index; [-]
R—correlation coefficient; [-]
τ—shear stress; [Pa]
τi- —shear stress measured at i-th rotational speed; [Pa]
Energies 2020, 13, 3192 17 of 17

τy —yield point; [Pa]

τ̂—shear stress determined from rheological model; [Pa]
τ—average value of shear stress; [Pa]
w/c—water cement ratio. [-]

References
1. API RP 13D. Rheology and Hydraulics of Oil-well Drilling Fluids; Norm of American Petroleum Institute:
Washington, DC, USA, 2006.
2. Bourgoyne, A.T.; Milheim, K.K.; Chenevert, M.E.; Young, F.S. Applied Drilling Engineering. In SPE Textbook;
SPE: Houston, TX, USA, 1986; 4.8.
3. Klotz, J.A.; Brigham, W.E. To Determine Herschel-Bulkley Coefficients. J. Pet. Technol. 1998, 50, 80–81.
[CrossRef]
4. Wiśniowski, R. Zastosowanie modelu Herschela-Bulkleya w hydraulice płuczek wiertniczych. Nowoczesne
Techniki i Technologie Bezwykopowe. Zeszyt nr 2/2000 Kraków 2000, 20–28.
5. Bair, S. Actual Eyring Models for Thixotropy and Shear-Thinning: Experimental Validation and Application
to EHD. J. Tribol. 2004, 126, 728–732. [CrossRef]
6. Pellegrino, C.; Faleshini, F. Sustainability Improvements in the Concrete Industry: Use of Recycled Materials for
Structure Concrete Production; 2016 Technology & Engineering; Springer International Publishing: Cham,
Switzerland, 2016; pp. 37–43.
7. Vom Berg, W. Influence of specific surface and concentration of solids upon the fow behaviour of cement
pastes. Mag. Concr. Res. 1979, 31, 211–216. [CrossRef]
8. Rawlings, J.O.; Pentula, S.G.; Dickey, D.A. Applied Regression Analysis: A Research Tool, 2nd ed.; Springer:
New York, NY, USA, 1998.
9. Azure, I.; Aloliga, G.; Doabil, L. Comparative Study of Numerical Methods for Solving Non-linear Equations
Using Manual Computation. Math. Lett. 2019, 5, 41–46. [CrossRef]
10. Bazaraa, M.S.; Sherali, H.S.; Shetty, C.M. Nonlinear Programing—Theory and Algorithms, 3rd ed.; John Wiley &
Sons, Inc.: Hoboken, NJ, USA, 2006.
11. Luenberger, D.G.; Yinyu, Y. Linear and Nonlinear Programming, 4th ed.; Springer Science+Business Media,
LLC: New York, NY, USA, 2008; pp. 14–19.
12. Boyd, S.; Vandenberghe, L. Convex Optimization; United States of America by Cambridge University Press:
New York, NY, USA, 2004.
13. Boyou, N.V.; Ismail, I.; Sulaiman, W.R.W.; Sharifi Haddad, A.; Husein, N.; Hui, H.T.; Nadaraja, K. Experimental
investigation of hole cleaning in directional drilling by using nano-enhanced water-based drilling fluids.
J. Pet. Sci. Eng. 2019, 176, 220–231. [CrossRef]
14. Sharma, M.M.; Zhang, R.; Chenevert, M.E. A new family of nanoparticle-based drilling fluids. In Proceedings
of the SPE Annual Technical Conference and Exhibition, San Antonio, TX, USA, 8–10 October 2012; pp. 1–13.
15. Smith, S.R.; Rafati, R.; Sharifi Haddad, A.; Cooper, A.; Hamidi, H. Application of aluminium oxide
nanoparticles to enhance rheological and filtration properties of water-based muds at HPHT conditions
Colloids Surfaces A Physicochem. Eng. Asp. 2018, 537, 361–371. [CrossRef]
16. Rafati, R.; Smith, S.R.; Sharifi Haddad, A.S.; Novara, R.; Hamidi, H. Effect of nanoparticles on the modifications
of drilling fluids properties: A review of recent advances Colloids Surfaces A Physicochem. Eng. Asp. 2018,
161, 61–76.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).