Advanced iterative procedures for solving the implicit Colebrook equation for fluid flow friction

Empirical Colebrook equation from 1939 is still accepted as an informal standard to calculate friction factor during the turbulent flow through pipes from smooth with almost negligible relative roughness to the very rough inner surface. The Colebrook equation contains flow friction factor in implicit logarithmic form where it is, aside of itself, a function of the Reynolds number Re and the relative roughness of inner pipe surface. To evaluate the error introduced by many available explicit approximations to the Colebrook equation, it is necessary to determinate value of the friction factor from the Colebrook equation as accurate as possible. The most accurate way to achieve that is using some kind of iterative methods. Usually classical approach also known as simple fixed point method requires up to 8 iterations to achieve the high level of accuracy, but does not require derivatives of the Colebrook function as here presented accelerated Householder approach (3rd order, 2nd order: Halley and Schroder method and 1st order: Newton-Raphson) which needs only 3 to 7 iteration and three-point iterative methods which needs only 1 to 4 iteration to achieve the same high level of accuracy. Strategies how to find derivatives of the Colebrook function in symbolic form, how to avoid use of the derivatives (Secant method) and how to choose optimal starting point for the iterative procedure are shown. Householder approach to the Colebrook equations expressed through the Lambert W-function is also analyzed. One approximation to the Colebrook equation based on the analysis from the paper with the error of no more than 0.0617% is shown.


Introduction
To evaluate flow resistance in turbulent flow through pipes, from smooth to rough, the empirical Colebrook equation is in common use (1) {Colebrook 1939}: (1) In the Colebrook equation λ represents Darcy flow friction factor, Re Reynolds number and ε/D relative roughness of inner pipe surfaces (all three quantities are dimensionless).
The experiment performed by Colebrook and White {Colebrook and White 1937} dealt with flow of air through one pipe (diameter D=53.5mm and length L=6m) with six different roughness of inner surface of the pipe artificially simulated with various mixtures of two sizes of sand grain (0.035mm and 0.35mm diameter) to simulate conditions of inner pipe surface from almost smooth to very rough. The sand grains were fixed using a sort of bituminous adhesive waterproof insulating compound to form five types of relatively uniform roughness of inner pipe surfaces while the sixth one was without sand, i.e. it was left smooth. The experiment revealed, contrary to the previous, that flow friction, λ does not have a sharp transition from the smooth to the fully rough law of turbulence. This evidence Colebrook {Colebrook 1939} later captured in the today famous and widely used empirical equation; Eq. (1).
The Colebrook function relates the unknown flow friction factor λ as function of itself, the Reynolds number Re and the relative roughness of inner pipe surface ε/D; λ=f(λ, Re, ε/D). It is valid for 4000<Re<10 8 and for 0<ε/D<0.05. The Colebrook equation is transcendent and thus cannot be solved in terms of elementary functions {Sonnad and Goudar 2007, Brkić 2011a, Mikata and Walczak 2016. Although empirical, and therefore with questionable accuracy its precise solution is sometimes essential in order to repeat or to evaluate previous findings in a concise way {Schockling et al. 2007, Allen et al. 2007, Langelandsvik et al. 2008.  Herwig et al. 2008, Brkić 2012a but which can be used as initial starting point for all cases covered by the Colebrook equation λ 0 =f(ε/D)→Eq.

2.b.) Householder's iterative methods.
On the other hand, Newton's method (also known as the Newton-Raphson method {Cajori 1911, Ypma 1995 needs only 3 to 7 iterations to reach the same level of accuracy. A shortcoming of Newton's method is that it additionally requires the first derivative of the Colebrook's function (here we show analytical form of the first derivative including the symbolic form generated in MATLAB {Shampine 2007}). Also knowing that the Newton-Raphson method is 1 st order of Householder's method {Householder 1970, Griffiths and Smith 2006}, we analyze 2 nd order which is known as Halley's {Halley 1694} and Schröder's {Schröder 1870 method, and also 3 rd order.
The 3 rd order methods use the third, the second and the first derivative, 2 nd order the second and the first while 1 st order use only the first derivative. Today, all mentioned types of iterative solutions can easily be implemented in software codes and they has been accepted as the most accurate way for solving the Colebrook equation and hence they are preferred compared to the graphical solution.
(2.c.) Three-point iterative methods. Three-point iterative methods needs only 1 iteration in three points x 0 , y 0 and z 0 (three internal iterations) to achieve high level of accuracy {Džunić et al. 2011, Petković et al. 2014, Sharma and Arora 2016 x 0 is initial starting point, y 0 is auxiliary step, while z 0 is the solution. Three-point methods are very accurate and can reach high accuracy some cases even after 1 to 2 iterations. Also slightly less accurate two-point methods exist.  Zigrang and Sylvester 1985, Giustolisi et al. 2011, Genić et al. 2011, Brkić 2011b, Brkić 2012b, Winning and Coole 2013. Numerous explicit approximations to the Colebrook equation are available in literature {Gregory and Fogarasi 1985, Zigrang and Sylvester 1985, Brkić 2011b, Genić et al. 2011, Brkić 2012b, Winning and Coole 2013 Iterative solution as the most accurate method is used for evaluation of accuracy of such approximations. Also, based on our findings, we provide an approximation; Eq. (28)  We also present distribution of the relative error in respect of the presented approximation over the applicability domain of the Colebrook equation.

Initial estimate of starting point for the iterative procedures
The starting point is a significant factor in convergence speed in three-point and Householder's method {Kornerup and Muller 2006} and there are different methods to choose a good start but here we examine 1) Starting point as function of input parameters, and 2) Initial starting point with the fixed value.
One of the essential issues in every iterative procedure is to choose good starting point {Moursund 1967, Taylor 1970}. Here we try to find the fixed starting point (initial value of the flow friction factor λ 0 or the related transmission factor ) valid for all cases from the practical domains of applicability of the Colebrook equation which is; Reynolds number Re; 4000<Re<10 8 and the relative roughness ε/D; 0<ε/D<0.05. In the cases when this approach does not work efficiently we show how to choose the starting value in function of the Reynolds number Re and the relative roughness ε/D, i.e. using some kind of rough approximation to the Colebrook equation which can be relatively inaccurate but simply and which put the initial value reasonable close to the final and accurate solution. This initial guess then need to be plugged into the shown numerical methods and iterated recursively few times (usually two or three times) to converge upon the final solution. In any case, a sample of size 65536 was considered for analysis of iteration methods. The input sample was generated according to the uniform density function of each input variable. The low-discrepancy Sobol sequences were employed {Sobol et al. 1992}. These so called quasirandom sequences have useful properties. In contrary to random numbers, quasirandom numbers covers the space more quickly and evenly. Thus, they leave very few holes.
The Colebrook equation also can be expressed in term of the Lambert W-function analytically; λ=f(λ, Re, ε/D)→λ=W(Re, ε/D) {Keady 1998, Goudar 2004, Brkić 2012c} (2) The initial starting point obtained using the previous equation is referred as "traditional", it introduces the maximal relative error of 80% over the domain of applicability of the Colebrook equation where the error can be neglected in case of fully developed turbulent flow through the pipes with very rough inner surface. To reach accuracy of λ i+1 -λ i ≤10 -8 usually 6 steps are enough regarding the Newton-Raphson method ( Figure 1).

Fixed initial starting point
An idea from geometry to find "center of gravity" is used for the points for which the Newton-Raphson, Hallaey's, Schröder's and three-point methods converge slowly. If we put the initial starting point in this zone, the less number of iterations is required to reach the final solution {Brkić 2014}.

Fixed initial starting point for the Newton-Raphson method
The "center of gravity" for the "slow area" in which the Newton-Raphson method requires increased number of iterations is shown in Figure 1. The "center of gravity" has coordinates: log(Re)=4.4322→Re≈27000 and -log(ε/D)=5.7311→ε/D≈1.85·10 -6 for which the flow friction factor λ and the corresponding transmission factor can be calculated using any of the available methods. In that way calculated x, became the starting point x 0 for all combinations of the Reynolds number Re and the relative roughness ε/D in the domain of applicability of the Colebrook equation. With this new starting point x 0 the maximal required number of iterations is 4 (Figure 3), while before in the worst case was 6 ( Figure 1) when the starting point x 0 was obtained through the "traditional formula" for this purpose; Eq.
(2), all valid for the case when the flow friction factor λ is calculated with accuracy of λ i+1λ i ≤10 -8 using Colebrook equation solved in Newton's procedure when calculation goes through the transmission factor x.
The physical interpretation of this "slow area" is in the fact that this area corresponds to the initial zone of turbulent flow through smooth pipes while Eq. (2) is accurate only for the fully developed turbulent flow through rough pipes. So, Eq. (2) can already obtain accurate solution in the case of fully developed turbulent flow through rough pipes even without the iterative process, where Eq. (2) introduces the relative error of almost 80% in the case of initial phases of turbulent flow through smooth pipes.
With the initial starting point fixed at the "centre of gravity" of the "slow area", in the worst cases, maximum 4 iterations as shown in Figure 3 are enough for the required accuracy of 10 -8 (before with the "traditional" version of initial value provided using Eq. (2) was 6 as indicated in Figure 1). The new fixed starting point is set as λ 0 =0.024069128765100981, i.e. x 0 =6.44569593948452. It corresponds to log(Re)=4.4322→Re≈27000 and -log(ε/D)=5.7311→ε/D≈1.85·10 -6 . The new starting point x 0 =6.44569593948452 is very robust and it seems to be an optimal starting point for all combinations when calculation goes through the transmission factor as explained in Section 3.2.

Fixed initial starting point for the Halley and Schröder method
Starting point for calculation through the Halley and Schröder method using Eq.
(2) requires in the worst cases up to 4 iterations to reach the required accuracy ( Figure 4). Compared with the Newton-Raphson method it is improvement of two iterations; up to 6 iterations required in Figure 1 and up to 4 iterations in Figure 4. The "worst-case" area for Halley and Schröder's method that requires 4 iterations using staring point Eq.

Starting point for the Lambert W-expressed Colebrook equation
The friction factor λ in the Colebrook equation can be expressed in explicit way through the Lambert Wfunction { Keady 1998, Sonnad and Goudar 2006, Sonnad and Goudar 2007, Clamond 2009, Brkić 2011a, Brkić 2012b. The Lambert W-function further can be evaluated using some types of Householders iterative methods as shown in Section 3.5 of this paper. The argument of the Lambert W-function in this case, exp(α 0 ) depends on both Re and ε/D, but as explained due to exponential form the calculation is not always possible and because of that limited possibility of use the appropriate starting point in this case is not evaluated {Goldberg 1991, Sonnad and Goudar 2004, Kornerup and Matula 2010.

Iterative Methods Adopted for the Colebrook Equation
Householder's method {Householder 1970} is a numerical algorithm for solving the nonlinear equation such as Colebrook's. During the Householder's procedure, in successive calculation, i.e. in iterative cycles the original assumed value of the unknown quantity (initial starting point {Kornerup and Muller 2006}) needs to be brought as much as possible close to the real value of the quantity using the least possible number of iterations. The same situation is with the three-point methods {Sharma and Arora 2016}.
The following types of the method are used in this paper: 1 st order Householder's method (Newton-Raphson {Cajori 1911, 2 nd order (Halley {Gander 1985, and 3 rd order, as well three-point methods {Sharma and Arora 2016}. All these methods require calculation of the derivatives which is usually underlined as the most important shortcoming of Householder's and three-point methods compared with the simple fixed point procedure in respect of the Colebrook equation. The Newton-Raphson and three-point method require only the first derivative, the Halley and Schröder the first and the second derivative, while the 3 rd requires the first, the second and the third derivative. In addition to the first derivate in analytical form, all required derivatives of the Colebrook function we present also in a simple and computationally inexpensive symbolic form. The derivatives in symbolic form were generated in MATLAB. In addition, the Secant method which does not require derivatives is shown as a

Direct calculation of λ with derivative calculated in analytical way
The proposed technique requires the Colebrook equation in the form f(λ, Re, ε/D)=0; Eq. (5) where λ is treated as variable, the first derivative f'(λ, Re, ε/D); Eq. (6) of the Colebrook equation with respect of λ and the initial value of the friction factor λ 0 as starting point. Most probably, the function will have residue f/f'≠0 which needs to be minimized through the iterative process.
Here are the required steps for the Newton-Raphson procedure: -The Colebrook equation in the form f(λ, Re, ε/D)=0; (5): -The first derivative f' with respect of λ in exact analytical way (6): -Initial value λ 0 selected as explained in Section 2 of this paper in order to calculate the residue f/f' and start the iterative procedure; Eq. (7): The procedure λ i+1 =λ i -f(λ i )/f'(λ i ) needs to be followed until the residue f(λ i )/f'(λ i )≈0.
(2), calculation of λ -Eq. (7), analytical derivative f'(λ) -Eq. (6) Re=5·10 6 , ε/D=2.5·10 -5 f(λ); Eq (5) f'(λ); Eq. (6)   Comparing the same approach but with the different starting points (Table 1 and 2), we can conclude that one single calculated negative value for flow friction factor λ can increase number of required iterations significantly. These negative values can occur if the initial starting point is chosen too far away from the final calculated solution. This problem can be overwhelmed with the Colebrook function slightly rearranged as in Section 3.2.

Indirect calculation of λ through the transmission factor
The appropriate form of the function is essential to reduce number of required iteration to reach the final solution. In order to accelerate the procedure an appropriate shift is used to provide some kind of linearization of the problem.
The Newton-Raphson procedure with these changes has similar steps as already shown: - The first derivative of Eq. (8) in respect to the transmission factor x can be calculated analytically, but also in symbolic form (where both approaches give identical results); Sections 3.2.1 and 3.2.2.

Indirect calculation of λ through the transmission factor with the derivative calculated
analytically -The first derivative f' with respect of x can be obtained analytically; Eq. (9), (also Eq. (11) gives the same results): -Initial value of the flow friction factor λ 0 should be chosen and the residue f/f' calculated in order to start the Newton-Raphson procedure; Eq. (10): The procedure x i+1 =x i -f(x i )/f'(x i ) should be followed until the residue f(x i )/f'(x i )≈0. Then the final solution is , where n=i+1 is the final iteration. Approach with the indirect calculation of λ through the transmission factor x is much more stable compared with the direct calculation of λ as can be seen from Tables 2 and 3 comparing the number of required iterations to reach the same accuracy (11 iterations for the direct approach compared with only 3 iterations in the indirect approach using fixed starting point x 0 =6.445695939 for Re=5·10 6 , ε/D=2.5·10 -5 ).

Indirect calculation of λ through the transmission factor with the symbolic derivative
-The exact analytical expression of the first derivative f' with respect of x can be obtained in MATLAB; Eq. (11), results are the same as using Eq. (9): -Initial value of the flow friction factor λ 0 should be chosen and the residue f/f' calculated in order to start the Newton-Raphson procedure; Eq. (12): The procedure x i+1 =x i -f(x i )/f'(x i ) should be followed until the residue f(x i )/f'(x i )≈0. Then the final solution is , where n=i+1 is the final iteration.
The iterative procedure can be accelerated using Halley's formula instead of the Newton-Raphson; Eq. (13): (13) In general x 1 =x i+1 and x 0 =x i ; i=0 to n, where n+1 is the final iteration in which x n ≈x n+1 .
The second derivative f''(x) in respect of x is required; Eq. (14): The Newton-Raphson method belongs to the 1 st order, and the Halley to the 2 nd order of Householder's method while the 3 rd order can be expressed using Eq. (15): (15) Again, x 1 =x i+1 and x 0 =x i ; i=0 to n, where n+1 is final iteration in which x n ≈x n+1 .
The required 3 rd derivative f'''(x) can be expressed using Eq. (16) Further x 1 =x i+1 and x 0 =x i ; i=0 to n, where n+1 is final iteration in which x n ≈x n+1 .

Secant method
Secant method is similar to the Newton-Raphson, it requires two starting points λ 0 and λ -1 but doesn't require calculation of the derivatives {Varona 2002}. The approach with the direct calculation of λ with the two required starting points λ 0 and λ -1 is formulated; Eq. (18): Counter i starts from i-1 and goes to n+1 in which λ n =λ n+1 .
The approach through the transmission factor x is; Eq. (19): As already described, counter i also starts from i-1 and goes to n+1 in which x n =x n+1 .
Flow friction factor λ is calculated in Tables 8 and 9 for two pairs of the Reynolds number and the relative roughness 1) Re=5·10 6 , ε/D=2.5·10 -5 and 2) Re=3·10 4 , ε/D=9·10 -3 using the Secant procedure with direct calculation of λ and indirect through the transmission factor x. Calculation using the Secant procedure, also confirms that the indirect calculation of λ through the transmission factor x requires in general less number of iterations to reach the same level of accuracy. In the case from Tables 8 and 9, the required number of iterations is 6 in direct calculation and 5 in indirect for Re=5·10 6 , ε/D=2.5·10 -5 , the required number of iterations is 4 in direct calculation and 3 in indirect for Re=3·10 4 , ε/D=9·10 -3 .

Expressed through the Lambert W-function
The Lambert W-function W(y) {Corless et al. 1996} is the solution of , which needs to be in the appropriate form; Eq.
Further in all cases, the Newton-Raphson, the Halley, and the Schröder; z 1 =z i+1 and z 0 =z i ; i=0 to n, where n+1 is final iteration in which z n ≈z n+1 .
Argument of the Lambert W-function y, in our case defined by Eq. (3), is y=Re/2.18. Therefore it does not depend on the relative roughness ε/D but only on the Reynolds number Re. In Table 10, z is calculated in the iterative procedure using the Newton-Raphson and Halley method for Re=5·10 6 and Re=3·10 4 where initial starting point is set as z 0 =15 as recommended in Section 2.3 of this paper (for z 0 <8.814, the Newton-Raphson procedure cannot start).

Approximations -Simplified equations for engineering practice
Using optimal fixed initial starting point for the Halley and the Schröder method as explained in Section 2.2.2, the first iteration of the procedures from Section 3.2.2, simplification using the fact that the first derivative of the Colebrook function is almost always near zero; f'≈0; and using acceleration through Eq.
The shown procedure is efficient and does not require extensive computing resources since the accuracy of more than 0.69% can be reached using only two logarithmic forms, while very high accuracy of 0.0617% using only three logarithmic forms (in all cases exponents are only whole numbers) {Clamond 2009, Giustolisi et al. 2011.

Conclusion
Today, a fast but reliable approximation of pipeline hydraulic is needed also reliability modelling of water and gas networks where a large number of network simulations of random component failures and their combinations must be automatically evaluated and statistically analyzed {Brkić 2009, Spiliotis and Tsakiris 2010, Brkić 2011e, Brkić 2011f, Brkić 2012e, Brkić 2016b. Accurate, fast and reliable estimation of flow friction factors are essential for evaluation of pressure drops and flows in large network of pipes, because, for example, compression station failure can be approximated on the transmission level by user-defined logic rules obtained from hydraulic software {Brkić and Tanasković 2008, Praks et al. 2015, Praks et al. 2017, Badami et al. 2018. Iterative solutions and approximations for calculation of flow friction factor are implemented in software packages which are in common use in everyday engineering practice {Brkić 2016b}. So in this paper we analyzed selected iterative procedures in order to solve the Colebrook equation {Chun and Neta 2017, Zhanlav et al. 2017} and we found that two or three iterations of Halley and the Schröder methods are suitable for the required accuracy needed for the engineering practice, when the fixed initial starting point of Section 2.2.2 is applied. On the other hand, using three-point iterative method with the same initial conditions, the required high accuracy can be reached after only one iteration but using three internal steps {Džunić et al. 2011, Petković et al. 2014 Moreover, to simplified calculation for engineering use we can recommend: 1. Knowing that the Colebrook equation is used in engineering practice only in the limited domain of the Reynolds number Re between 4000 and 10 8 , and of the relative roughness of inner pipe surface ε/D up to 0.05, we evaluated number of iteration from that domain to reach sufficient accuracy and we detected zones in which additional number of iterations are required. Therefore we put the fixed initial point to start the iterative procedure in that zones in order to decrease number of required iterations. 2. Using the simplified Halley and the simplified Schröder procedure with the fixed starting point, after only one iteration we developed approximation that is with the error up to 8.29%. This is near the accurate value and therefore for the second iteration the simplified Newton-Raphson method can be used to reach accuracy of 0.69% and for extreme accuracy the third iteration need to be used to reach accuracy of 0.0617%. The simplification is in fact that the first derivative of the Colebrook function in our case is always near zero; f'≈0 (for f'→0 the Newton-Raphson method become the fixed-point method). Accuracy of 0.69% can be reached using only two logarithmic forms and of 0.0617% with only three logarithmic forms which does not require extensive computational efforts (the goal is to use the least possible number of logarithmic functions or functions with non-integer power) {Giustolisi et al 2011, Vatankhah 2018}. Moreover, the computational cost of iterations can be also reduced using Padé polynomials {Praks and Brkić 2018}.
We analyzed also the Colebrook equation expressed through the Lambert W-function and we found that the Halley and the Schröder methods can be advised in comparison to the Newton-Rapson method (where the problem with the initial starting point exists).

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper. The views expressed are those of the authors and may not in any circumstances be regarded as stating an official position of the affiliated employers of the authors; European Commission and Technical University Ostrava. Both authors contributed equally to this study.

Data availability statement
All conclusions are based on the calculation shown in this paper. All data which is required to repeat the calculation can be obtained using equation shown in this paper. Examples are also shown in Tables to support conclusions.