*Communication* **Rational Approximation for Solving an Implicitly Given Colebrook Flow Friction Equation**

**Pavel Praks 1,2,\* and Dejan Brki´c 1,2,3,\***


Received: 21 November 2019; Accepted: 13 December 2019; Published: 20 December 2019 -

**Abstract:** The empirical logarithmic Colebrook equation for hydraulic resistance in pipes implicitly considers the unknown flow friction factor. Its explicit approximations, used to avoid iterative computations, should be accurate but also computationally efficient. We present a rational approximate procedure that completely avoids the use of transcendental functions, such as logarithm or non-integer power, which require execution of the additional number of floating-point operations in computer processor units. Instead of these, we use only rational expressions that are executed directly in the processor unit. The rational approximation was found using a combination of a Padé approximant and artificial intelligence (symbolic regression). Numerical experiments in Matlab using 2 million quasi-Monte Carlo samples indicate that the relative error of this new rational approximation does not exceed 0.866%. Moreover, these numerical experiments show that the novel rational approximation is approximately two times faster than the exact solution given by the Wright omega function.

**Keywords:** hydraulic resistance; pipe flow friction; Colebrook equation; Colebrook–White experiment; floating-point computations; approximations; Padé polynomials; symbolic regression

#### **1. Introduction**

The Colebrook equation [1] for turbulent flow friction is implicitly given with respect to the unknown Darcy flow friction λ, as shown in Equation (1):

$$\frac{1}{\sqrt{\lambda}} = -2 \cdot \log\_{10} \left( \frac{2.51}{R\varepsilon} \cdot \frac{1}{\sqrt{\lambda}} + \frac{\varepsilon}{3.71} \right) \tag{1}$$

where:

λ—Darcy flow friction factor (dimensionless) *Re*—Reynolds number, 4000 < Re < 10<sup>8</sup> (dimensionless) ε—relative roughness of inner pipe surface, 0 < ε < 0.05 (dimensionless)

As a PhD student at the Imperial College in London, Colebrook developed his empirical equation based on the data from his joint experiment with his supervisor, Prof. White [2]. They experimented with flow of air through pipes with different roughness of the inner pipe surface. The experiment by Colebrook and White was described in a scientific journal and published in 1937 [2], while the related empirical equation by Colebrook was published in 1939 [1].

Compared with some other experimental findings [3], the Colebrook equation fits the friction factor within a few dozen percent of error [4]. The Colebrook equation over the last 80 years has been seen by the industry as an informal standard for flow friction calculation and has been very well accepted in everyday engineering practice. Based on the Colebrook equation, Moody developed a diagram that was used before the era of computers for graphical determination of turbulent flow friction [5]. Today, such nomograms have been replaced by explicit approximations, which introduce some value of error [6,7], or by iterative methods [8–10].

In a central computer processor (CPU), transcendental functions such as logarithmic, exponential, or functions with non-integer terms require execution of numerous floating-point operations, and therefore they should be avoided whenever possible [11–18]. Praks and Brki´c [19] recently developed a one-log call iterative method, which uses only one computationally demanding function, and even then only in the first iteration (for all succeeding iterations, cheap Padé approximants are used [20,21]). Based on that approach, few very accurate and efficient explicit approximations suitable for coding and for engineering practice have been constructed [22]. In addition, the same authors developed few approximations of the Colebrook equation based on the Wright ω-function, which are among the most accurate to date [23–25]. On the other hand, they contain one or two logarithmic functions, depending on the chosen version [23,24] (these procedures are based on the previous efforts by Praks and Brki´c for symbolic regression [26] and by Brki´c with Lambert W-function [27–29]).

In this communication, we make a step forward and we offer for the first time a procedure for the approximate solution of the Colebrook equation based only on rational functions. The presented novel rational approximation procedure introduces a relative error of no more than 0.866% for 0 < ε < 0.05 and 4000 < Re < 10<sup>8</sup> (as used in engineering practice). The rational approximation procedure is suitable for computer codes (open-source code in commercial Matlab 2019a is given in this communication, and in addition, it is compatible with freeware GNU Octave, version 5.1.0).

After introductory Section 1, Section 2 of this communication gives a short overview of mathematical methods used for the proposed rational approximation procedure. Section 3 describes the rational approximation procedure in detail (including error analysis), Section 4 provides software code along with the algorithm to be followed for the rational approximation approach, while Section 5 contains concluding remarks.

#### **2. Mathematics Behind the Proposed Approximation**

Our rational approximation approach is based on Padé approximants [20,21], symbolic regression [30,31] and although not used directly, it is inspired by the Wright ω-function, a cognate of the Lambert W-function [32]. To avoid detailed explanations about the Lambert W-function [33], here it should be noted that in this context it is used to transform the Colebrook equation from the shape implicitly given in respect to the unknown flow friction factor to the explicit form [23,24,27–29,34,35].

#### *2.1. Padé Approximants*

The ratio of two power series with properly chosen coefficients of the numerator and denominator can approximate very accurately various functions in a narrow zone around the chosen expanding point. For the expressions in the numerator and denominator, Padé approximants [20,21] use rational functions of given order instead of serial expansions. So, in other words, the Padé approximants can estimate functions usually in a narrow zone as the quotient of two polynomials, often has better approximation properties compared with its truncated Taylor series. Being a quotient, the Padé approximants are composed of lower-degree polynomials, where the degree of polynomials can be chosen according to needs. We use Matlab 2019a in order to generate the needed Padé approximants as replacements of the logarithmic function in our rational approximation approach. We do not use expressions with non-integer exponents, because according to Clamond [10], in the software interpretation it is evaluated through one exponential and one logarithmic function (for example *B*<sup>κ</sup> = *e*κ· ln (*B*), where κ is in most cases a non-integer). The computational complexity of an algorithm describes the amount of resources required to run it, for example the execution time. Winning and Coole [13] performed 100 million calculations for each mathematical operation using random inputs, with each repeated five times, and found that the most efficient operation for addition requires 23.4 s. According to them, relative effort for computation referred to addition-1 as a reference for the following values: logarithm to base, 10–3.37; fractional exponential, 3.32, cubed root, 2.71; natural logarithm, 2.69; cubed, 2.38; square root, 2.29; squared, 2.18; multiplication, 1.55, division, 1.35; subtraction, 1.18.

#### *2.2. Symbolic Regression*

Symbolic regression is a machine learning approach for finding approximate functions based on evolutionary or genetic algorithms [36]. To avoid imposing prior assumptions on the model, symbolic regression has the ability to search through the space of mathematical expressions to look for an approximate function that best fits a given dataset [31].

We used Eureqa [30], a symbolic regression engine, to obtain our final model. HeuristicLab [37], a software environment for heuristic and evolutionary algorithms, including symbolic regression, can be used instead.

#### **3. Routine Based on Polynomial-Form Expressions**

#### *3.1. Replacement of Logarithmic Function*

The Colebrook equation can be accurately approximated using a rational approximation procedure, as shown in Equation (2):

$$\frac{1}{\sqrt{\lambda}} \approx -0.8686 \cdot (\zeta\_1 + \zeta\_2) \tag{2}$$

where:

$$\zeta\_1 = 0.02087 \cdot r - 0.07659 \cdot p(r) - \frac{0.5994}{p(r) + 3.846} - \frac{0.0007232}{r} - 0.00007489 \cdot r^2 + 0.1391$$

$$\zeta\_2 = p(r) - 7.93$$

$$p(r) = \frac{r \cdot (r \cdot (11 \cdot r + 27) - 27) - 11}{r \cdot (r \cdot (3 \cdot r + 27) + 27) + 3}$$

$$r = 2777.77 \cdot \left(\frac{2.51 \cdot p\_0}{Re} + \frac{\varepsilon}{3.71}\right)$$

$$p\_0 = \frac{2600 \cdot Re}{657.7 \cdot Re + 214600 \cdot Re \cdot \varepsilon + 12970000} - 13.58 \cdot \varepsilon + \frac{0.0001165 \cdot Re}{0.00002536 \cdot Re \cdot \varepsilon + Re \cdot \varepsilon + 105.5} + 4.227$$

Consequently, Equation (2) contains only rational functions, where:

<sup>ζ</sup><sup>1</sup> <sup>+</sup> <sup>ζ</sup><sup>2</sup> rational approximation of ln 2.51 *Re* · <sup>1</sup> √ <sup>λ</sup> <sup>+</sup> <sup>ε</sup> 3.71 ;

ζ<sup>1</sup> a rational function that corrects error caused by Padé approximant *p*(*r*);

<sup>ζ</sup><sup>2</sup> shifted Pad<sup>é</sup> approximant *<sup>p</sup>*(*r*), where the shift <sup>−</sup>7.93 <sup>≈</sup> *ln*(0.00036) where 0.00036 <sup>≈</sup> <sup>1</sup> 2777.77 ;


and where: <sup>−</sup><sup>2</sup> ln(10) <sup>≈</sup> <sup>−</sup><sup>2</sup> 2.302585093 ≈ −0.868, as log10(ς) <sup>=</sup> ln(ς) ln(10).

Function <sup>ζ</sup><sup>2</sup> approximates the required logarithmic function ln 2.51 *Re* · <sup>1</sup> √ <sup>λ</sup> <sup>+</sup> <sup>ε</sup> 3.71 , while ζ<sup>1</sup> corrects its error, as *r* is not always close to the expansion point, because of the large variability of input parameters of the Colebrook equation (Equation (1)). The rational function ζ<sup>1</sup> and also the starting point *p*<sup>0</sup> were found by symbolic regression software Eureqa [31], whereas the shift −7.93 in ζ<sup>2</sup> was found in order to minimize the error of the Padé approximation of *p*(*r*) ≈ ln(*r*) for the Colebrook equation, as ln(0.00036·*r*) <sup>≈</sup> ln(*r*) <sup>−</sup> 7.93, where 0.00036 <sup>≈</sup> <sup>1</sup> 2777.77 . Variable precision arithmetic (VPA) at 4 decimal digit accuracy is assumed for ζ<sup>1</sup> and for *p*0. The Padé approximant *p*(*r*) is given in

Horner nested polynomial form generated in Matlab 2019a. It is of order /2,3/, which means that the polynomial in numerator contains a monomial of highest degree 2, while the denominator of degree 3. Any other suitable Padé polynomial of any other degree in any other form that can substitute the natural logarithm around the needed expansion point can be used, but in such cases, the rational approximation procedure should be tested again, because these changes can affect the value of the final relative error and its distribution.

#### *3.2. Error Analysis*

Distribution of the relative error for the proposed rational approximation procedure, Equation (2), is given in Figure 1. The maximal relative error goes up to 0.866% for 0 < ε < 0.05 and 4000 < Re < 10<sup>8</sup> (as used in engineering practice). The highest error using 2 million input pairs found by the Sobol quasi-Monte Carlo method is for Re <sup>=</sup> 71987 and <sup>ε</sup> <sup>=</sup> 3.1711·10−<sup>7</sup> [38]. The Colebrook equation is empirical and it follows logarithmic law, so for the procedure with only rational functions, this level of error is acceptable [39]. For example, Pimenta et al. [7] classify the approximation of Sonnad and Goudar [40] with relative error of up to 3.17%. With two logarithms and one non-integer power, this method [40] belongs to the group of approximations with higher performance indexes and precision. Brki´c [6] estimates the relative error of Sonnad and Goudar [40] to be up to 0.8%, which is similar error compared with the rational approximation approach presented here; the same methodology as in Brki´c [6] is used for Figure 1.

**Figure 1.** The distribution of the relative error for the proposed rational approximation.

The error in the proposed rational approximation can be possibly reduced by optimizing numerical values of parameters using a genetic algorithms [41] approach with the methodology described by Brki´c and Cojbaši´ ´ c [42]. However, the presented rational approximation approach has too many numerical parameters, meaning such an optimization would be very complex. Further simplifications rather should go in the direction of simplification of *p*0, ζ<sup>1</sup> and ζ2, but keeping the same or increasing accuracy.

#### *3.3. Computational Costs*

The efficiency of the proposed procedure is tested using 2 million input pairs found by the Sobol quasi-Monte Carlo method [38]. The tests were performed using Matlab R2019a. The tests revealed that Equation (2) needs 0.56 s to calculate the friction factor λ for 2 million input pairs, or for the Reynolds number *Re* and the roughness of the inner pipe surface ε. On the other hand, the exact solution given by the Wright ω-function [23] implemented by the Matlab library "wrightOmegaq" [43] took 1.1 s for the same 2 million the tested pairs.

Our novel rational approximation of the Colebrook equation given with Equation (2) is approximately two times faster than the exact approach using the Wright ω-function [23], as the speed ratio is 1.1/0.56–1.96. On the other hand, the approach with the Wright ω-function [23] gives the exact solution, which requires two logarithms and one Wright ω-function, while this communication presents a rational approximation.

#### **4. Software Description**

The presented rational approximation approach for solving Colebrook's equation for flow friction was thoroughly tested at IT4Innovations, National Supercomputing Center, VŠB-Technical University of Ostrava, Czech Republic.

#### *4.1. Algorithm*

The simple algorithm of the rational approximation of the Colebrook equation is presented in Figure 2. The algorithm contains only one branch and is without loops.

**Figure 2.** Algorithm of the proposed rational approximation of the Colebrook equation.

#### *4.2. Open-Source Software Code*

The code is given in Matlab format, which is compatible with the freeware GNU Octave, but it can be easily transposed in any programming language (input parameters: R is the Reynolds number Re, K is the relative roughness of inner pipe surface ε; output parameter; L is the Darcy flow friction factor λ):

*x* = (2600\*R)/(657.7\*R + 214600\*R.\*K + 1.297e+7)−13.58\*K + (1.165e-4\*R)/(2.536e-5\*R + R.\*K + 105.5) + 4.227 *y*<sup>0</sup> = 2.51\*x./R + K./3.71; r = y0\*2777.77 *p*<sup>r</sup> = (r.\*(r.\*(11\*r + 27) − 27) − 11)/(r.\*(r.\*(3\*r + 27) + 27) + 3) *k*<sup>1</sup> = @(r) 0.02087\*r− 0.07659\*pr− 0.5994/(pr + 3.846) − 7.232e-4./r−7.489e-5\*r.ˆ2 + 0.1391 *x* = −0.8686\*(k1(r) + pr − 7.93)

$$L = 1/\infty. \text{ ?}$$

#### **5. Conclusions**

We provide a novel rational approximation procedure for solving the logarithmic Colebrook equation for flow friction. Instead of transcendental functions (logarithms, non-integer power) that are used in the classical approach, in this communication, we replace the logarithm with its Padé approximant and a simple rational function, which was found using artificial intelligence (symbolic regression), in order to minimize the error. Although the new rational approximation may seem unintelligible to human eyes, results of 2 million input pairs found by the quasi-Monte Carlo method [38] confirm that the relative error of this new approximation does not exceed 0.866%, which is acceptable for the empirical Colebrook law [44] (trade-off between model complexity and accuracy [45,46]). Consequently, numerical experiments on 2 million of quasi-Monte Carlo pairs indicates that the rational approximation presented here provides for Colebrook's flow friction model a useful combination of Padé approximants and artificial intelligence (symbolic regression).

**Author Contributions:** P.P. got the idea for the presented rational approximation approach and developed its first version. D.B. put the rational approximation approach into a form suitable for everyday engineering use. D.B. made a draft of this short note and prepared figures. P.P. wrote the software code in consultations with D.B. and tested it. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work of P.P. has been partially supported by the Technology Agency of the Czech Republic under the project "National Centre for Energy" TN01000007 and "Energy System for Grids" TK02030039. The work of D.B. has been partially supported by the Ministry of Education, Youth, and Sports of the Czech Republic through the National Programme of Sustainability II (NPS II), IT4Innovations Excellence in Science project LQ1602, and by the Ministry of Education, Science, and Technological Development of the Republic of Serbia under the project "Development of new information and communication technologies, based on advanced mathematical methods, with applications in medicine, telecommunications, power systems, protection of national heritage, and education" iii44006.

**Acknowledgments:** We especially thank Marcelo Masera for his scientific supervision and approval. We acknowledge support from the European Commission, Joint Research Centre (JRC), Directorate C: Energy, Transport, and Climate, Unit C3: Energy Security, Distribution, and Markets.

**Conflicts of Interest:** The authors declare no conflict of interest. Neither the European Commission, the VŠB, Technical University of Ostrava, the research and development Center "Alfatec", nor any person acting on their behalf is responsible for any use of this publication.

#### **Notations**

The following symbols are used in this Communication:

λ Darcy, Darcy–Weisbach, Moody, or Colebrook flow friction factor (dimensionless)


ω Wright ω-function (Wright omega function)

#### **References**


© 2019 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/).
