# Advanced Iterative Procedures for Solving the Implicit Colebrook Equation for Fluid Flow Friction

Authors: Pavel Praks; Dejan Brkić
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

The empirical Colebrook equation from 1939 is still accepted as an informal standard way to calculate the friction factor of turbulent flows (4000 < Re < 108) through pipes with roughness between negligible relative roughness (ε/D⟶0) to very rough (up to ε/D = 0.05). The Colebrook equation includes the flow friction factor λ in an implicit logarithmic form, λ being a function of the Reynolds number Re and the relative roughness of inner pipe surface ε/D: λ = f(λ, Re, ε/D). To evaluate the error introduced by the many available explicit approximations to the Colebrook equation, λ ≈ f(Re, ε/D), it is necessary to determinate the value of the friction factor λ from the Colebrook equation as accurately as possible. The most accurate way to achieve that is by using some kind of the iterative method. The most used iterative approach is the simple fixed-point method, which requires up to 10 iterations to achieve a good level of accuracy. The simple fixed-point method does not require derivatives of the Colebrook function, while the most of the other presented methods in this paper do require. The methods based on the accelerated Householder’s approach (3rd order, 2nd order: Halley’s and Schröder’s method, and 1st order: Newton–Raphson) require few iterations less, while the three-point iterative methods require only 1 to 3 iterations to achieve the same level of accuracy. The paper also discusses strategies for finding the derivatives of the Colebrook function in symbolic form, for avoiding the use of the derivatives (secant method), and for choosing an optimal starting point for the iterative procedure. The Householder approach to the Colebrook’ equations expressed through the Lambert W-function is also analyzed. Finally, it is presented one approximation to the Colebrook equation with an error of no more than 0.0617%.


---

## Body

## 1. Introduction

To evaluate flow resistance in turbulent flow through rough or smooth pipes, the empirical Colebrook equation is in common use [1]:(1)1λ=−2⋅log102.51/Re⋅λ+ε/3.7⋅D︸Colebrook.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 [2] dealt with flow of air through a pipe, diameter D = 53.5 mm, and length L = 6 m, with six different roughness of inner surface of the pipe artificially simulated with various mixtures of two sizes of sand grain (0.035 mm and 0.35 mm 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, that is, it was left smooth. The experiment revealed, contrary to the previous, that the flow friction, λ, does not have a sharp transition from the smooth to the fully rough law of turbulence. This evidence Colebrook [1] later captured in today famous and widely used empirical equation, Equation (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 < 108 and for 0 < ε/D < 0.05. The Colebrook equation is transcendent and thus cannot be solved in terms of elementary functions [3–6]. Although empirical, and therefore with questionable accuracy, its precise solution is sometimes essential in order to repeat or to evaluate the previous findings in a concise way [7–9].Few approaches are available today for solving the Colebrook equation:(1)
Graphical solution—Moody diagram: To represent the Colebrook equation graphically, Rouse in 1942 had developed an appropriate diagram which Moody later adapted in 1944 in the famous diagram widely used in the past in engineering practice [10, 11]. The diagram was preferred because the Colebrook equation is implicitly given. Today, graphical solution has only value for educational purposes.(2)
Iterative solution of the Colebrook equation:(a)
Simple fixed-point iterative method. The simple fixed-point iterative method [12] is in common use for solving accurately the Colebrook equation (special case of the Colebrook equation for Re ⟶ ∞ gives explicit form valid only for the fully turbulent flow in rough pipes [13–16] but which can be used as initial starting point for all cases covered by the Colebrook equation λ0 = f(ε/D)⟶Equation (2); now using the Colebrook equation, new value can be calculated λ1 = f(λ0; Re; ε/D); starting from i = 1, the procedure λi+1 = f(λi; Re; ε/D); i = i + 1 goes until λi ≈ λi+1, where we set λi+1 − λi ≤ 10−8). It usually reaches the satisfied accuracy after no more than 10 iterations.(b)
Householder’s iterative methods. On the other hand, Newton’s method (also known as the Newton–Raphson method [17–19]) needs few iterations less compared to the fixed-point method to reach the same level of accuracy. A shortcoming of Newton’s method is that it additionally requires the first derivative of the Colebrook function (here we show analytical form of the first derivative including the symbolic form generated in MATLAB [20]). Also knowing that the Newton–Raphson method is the 1st order of Householder’s method [21, 22], here we also analyze the 2nd order, which is known as the Halley [23] and the Schröder [24, 25] method, and also the 3rd order. The third-order methods use the third, the second, and the first derivative, the 2nd order use the second and the first, while the 1st order use only the first derivative. Today, all mentioned types of iterative solutions can easily be implemented in software codes and they are accepted as the most accurate way for solving the Colebrook equation, and hence, they are preferred compared to the graphical solution.(c)
Three-point iterative methods. Three-point iterative methods need only 1 to 3 iterations in three points x0, y0, and z0 (three internal iterations) to achieve the high level of accuracy [26–28]. x0 is initial starting point, y0 is auxiliary step, while z0 is the solution. Three-point methods are very accurate and can reach high accuracy in some cases even after 1 to 2 iterations. Also slightly less fast two-point methods in terms of required number of iterations to reach the demanded accuracy do exist.(3)
Approximations of the Colebrook equation: Colebrook’s equation can be expressed in the explicit form only in an approximate way: λ≈f(Re, ε/D) [29–36]. Numerous explicit approximations to the Colebrook equation are available in the literature [29–32, 34–36]. The iterative solutions as the most accurate methods are used for evaluation of accuracy of such approximations. Also, based on our findings, we provide an approximation, Equation (28), with the error of no more than 0.69% and 0.0617%. The Colebrook equation can also be approximately simulated using Artificial Neural Networks [37–39].(4)
Lambert W-function: Until now, the only one known way to express the Colebrook equation exactly in explicit way is through the Lambert W-function, λ = W(Re, ε/D) [3, 8, 40–43], where further evaluation of the Lambert W-function can be only approximated [44–48]. Here, we show the procedure how to solve the Lambert W-function using the Householder iterative procedure (2nd order: Halley’s method and 1st order: Newton–Raphson). Also approach with the shifted Lambert W-function in terms of the Wright Ω-function exists [40, 43].In this paper, we show the three-point and the Householder iterative procedures (the 3rd order, the 2nd order: Halley’s [49] and Schröder’s method, and the 1st order: Newton–Raphson) with the additional recommendations in order to solve the empirical Colebrook equation implicitly given in respect of the flow friction factor λ. The goal of this paper is to show the improved iterative solutions which can obtain the value of the unknown friction factor λ accurately after the least possible number of iterations. Additionally, we developed a strategy how to choose the best starting point [50] for the iterative procedure in the domain of interest of the Colebrook equation, how to generate required symbolic derivatives to the Colebrook equation in MATLAB, and how to avoid use of the derivatives (secant method). Finally, we use findings from our paper to present a novel explicit approximation of the Colebrook equation, which would be interesting for engineering practice. We also present distribution of the relative error in respect of the presented approximation over the applicability domain of the Colebrook equation.To evaluate the efficiency of the presented methods, the unknown flow friction factorλ is calculated for two pairs of the Reynolds number Re and relative roughness of inner pipe surfaces ε/D: (1) (Re = 5·106, ε/D = 2.5·10−5)⟶λ = 0.010279663295529 and (2) (Re = 3·104, ε/D = 9·10−3)⟶λ = 0.038630738574792.


## 2. Initial Estimate of Starting Point for the Iterative Procedures

The starting point is a significant factor in convergence speed in the three-point and the Householder methods [50], and there are the different methods to choose a good start, but here we examine (1) starting point as function of the input parameters and (2) initial starting point with the fixed value.One of the essential issues in every iterative procedure is to choose the good starting point [51, 52]. Here, we try to find the fixed starting point (the initial value of the flow friction factor λ0 or the related transmission factor x0=1/λ0) valid for all cases from the practical domains of applicability of the Colebrook equation which is for the Reynolds number Re, 4000 < Re < 108, and for 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, that is, using some kind of the rough approximations 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 needs to be plugged into the shown numerical methods and iterated recursively few times (usually two or three times and up to ten in the worst case) to converge upon the final solution. In any case, a sample of size 65536 was considered for analysis of the iteration methods. The input sample was generated according to the uniform density function of each input variable. The low-discrepancy Sobol sequences were employed [53]. These so-called quasirandom sequences have useful properties. In contrary to the random numbers, quasirandom numbers cover the space more quickly and evenly. Thus, they leave very few holes.The Colebrook equation can also be expressed in terms of the Lambert W-function analytically,λ = f(λ, Re, ε/D)⟶λ = W(Re, ε/D) [41, 42, 54]. The Lambert W-function further can be evaluated only approximately through the Householder iterative procedures which also require the appropriate initial starting point. The analysis of this initial starting point has wider applicability, because the Lambert W-function has extensive use in many branches of physics and technology [55, 56].


### 2.1. Starting Point as Function of Input Parameters

#### 2.1.1. Starting Point as Function of the Relative Roughnessε/D (When Re⟶∞)

The special case of the Colebrook equation when Re⟶∞ physically means that the flow friction factor λ in that case depends only on ε/D, for Re⟶∞, λ = f(ε/D), that is, the flow friction factor λ is not implicitly given [14]. In that way, the starting point can be calculated using the explicit equation which has only one variable, λ0 = f(ε/D) (Equation (2)). The results obtained in that way are accurate only for the case Re⟶∞ but for the smaller values of the Reynolds number Re which corresponds to the smooth turbulent flow, the error can goes up to 80% [13, 57]. Anyway, in that way, calculated value can be efficiently used as an initial starting guess for the iterative procedures for the whole domain of applicability of the Colebrook equation.(2)x0=1λ0=−2⋅log10ε/3.7⋅D︸roughpartofColebrook.The initial starting point obtained using the previous equation is referred as “traditional,” and 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 the accuracy ofλi+1 − λi ≤ 10−8, usually 6 steps are enough regarding the Newton–Raphson method (Figure 1).Figure 1
“Slow area” which requires 6 iterations to reach accuracy of 10−8 regarding the “traditional” option for the starting point calculated through Equation (2), for solving Colebrook’s equation in Newton’s procedure when calculation goes through the transmission factor x.


#### 2.1.2. Starting Point Obtained Using Approximations to the Colebrook Equation

Every approximation to the Colebrook equation,λ ≈ f(Re, ε/D), can be used to put the initial starting point as close as possible near the final accurate solution [30]. For example, using one of the approximations with the error of up to 10% for calculation of the initial starting point x0=1/λ, the final accurate value of the flow friction factor λ is reached in the worst-case scenario after 3 iterations using the Colebrook equation and one of the procedures from Section 3.2. After 3 iterations, the whole practical domain of applicability of the Colebrook equation is covered with the difference between the two final iterations less than 10−8, λi+1 − λi ≤ 10−8 (Figure 2). In average, the method requires 2.7 iterations in average for all cases with the set precision (stopping criterion) very close to zero (about 10−8) when calculation goes through the transmission factor x. The results from Figure 1 are from the 65536 pairs of the Reynolds number Re and the relative roughness ε/D over the domain of applicability of the Colebrook equation domain (values of the Reynolds number Re between 4000 and 108 and the relative roughness ε/D between 0 and 0.05, dividing them into 256 points each).Figure 2
Area in which 2 iterations (left) and 3 iterations (right) are sufficient to calculate the flow friction factorλ with accuracy of 10−8 using Colebrook equation solved in Newton’s procedure when calculation goes through the transmission factor x and using approximation with error of up to 10%.


### 2.2. Fixed Initial Starting Point

An idea from geometry to find “center of gravity” is used for the points for which the Newton–Raphson, the Halley, Schröder, and the 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 [58].


#### 2.2.1. Fixed Initial Starting Point for the Newton–Raphson Method

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


#### 2.2.2. Fixed Initial Starting Point for the Halley and Schröder Method

The starting point for calculation through the Halley and the Schröder method using Equation (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’s and Schröder’s method that requires 4 iterations using staring point Equation (2) has coordinates: (log10 (Re) = 5.3108⟶Re ≈ 204550; –log10 (ε/D) = 4.9431⟶ε/D ≈ 1.14·10−5)⟶λ0 = 0.015663210285978339, that is, x0 = 7.990256504. This value is the new optimal initial starting point in the case of the Halley and the Schröder method.Figure 4
“Slow area” which requires 4 iterations to reach accuracy of 10−8 regarding the “traditional” option for the starting point calculated through Equation (2), for solving Colebrook’s equation in Halley’s and Schröder’s procedure when calculation goes through the transmission factor x.With the new initial starting pointx0 = 7.990256504, three iterations are required at maximum to reach the required accuracy in case of the Halley and the Schröder method (Figure 5) as described in Section 3.2.Figure 5
Decreased maximal number of required iterations from 4 to 3 to reach accuracy of 10−8 for solving Colebrook’s equation in Halley’s and Schröder’s procedure when calculation goes through the transmission factor x where the initial starting point is with the fixed value: x0 = 7.990256504.


#### 2.2.3. Fixed Initial Starting Point for the Three-Point Iterative Methods

The optimal normalized parameters for the fixed initial starting point for the three-point iterative methods explained in Section3.4. of this paper [26–28] are as follows: (log10 (Re) = 4.90060379974617⟶Re = 79543.33576; –log (ε/D)  =  5.33355157079189 ⟶ε/D  =  4.63926·10−6)⟶λ0 = 0.018904186734624⟶x0 = 7.273124147. The Džunić–Petković method is shown in Section 3.4. Additional recommendations about the initial starting point regarding the three-point iterative methods can be found in [59].


### 2.3. Starting Point for the Lambert W-Expressed Colebrook Equation

The friction factorλ in the Colebrook equation can be expressed in the explicit way through the Lambert W-function [3, 5, 8, 30, 42, 60, 61]. The Lambert W-function can further be evaluated using some types of the Householder iterative methods as shown in Section 3.5 of this paper.The Colebrook equation in a closed form through the Lambert W-function can be expressed in two ways, Equations (3) and (4). The first expression is [3, 30, 60, 62] as follows:(3)1λ=−2⋅log102⋅2.51⋅WyRe⋅ln10+ε3.7⋅D=−2⋅log1010−Wy/ln10+ε3.7⋅D,where y=Re⋅ln10/2⋅2.51≈Re/2.18.The argument of the Lambert W-function in this case depends only on the Reynolds number,λ0 = Re/2.18. Knowing that the practical range of the Reynolds number goes from 4000 to 108, the argument of the Lambert W-function is goes from about 1835 to 45871560, where W(1835) = 5.763291081 and W(45871560) = 14.93748223. The Halley procedure is fast and any initial starting point can be chosen between 5 and 15, but the Newton–Raphson method is very slow and we found that, for the best results, the initial starting point 15 has to be chosen. Note that for W(45871560) = 14.93748223, the Newton–Raphson procedure does not work in Excel for the values of initial starting point lower than 8.814.Due to transformations of coefficients, Equation (3) can introduce the relative error up to 2% and should be considered as explicit approximation to the Colebrook equation rather to its equivalent [63].The second expression is as follows [42, 54]:(4)λ=2ln10⋅Weα−Re⋅ε/D3.7⋅2.51−2,where α=Re⋅ε/D⋅ln10/2⋅2.51⋅3.7−ln2⋅3.7/Re⋅ln10.Argument of the Lambert W-function in this case iseα which for the certain combinations of the Reynolds number Re and the relative roughness ε/D from the practical domain of the Colebrook equation is too big to be calculated in registers of computers [6, 41, 54]. This can be overwhelmed with the Wright Ω-function, ωβ=Weα [40, 43, 64–66].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 [41, 54, 67, 68].


## 2.1. Starting Point as Function of Input Parameters

### 2.1.1. Starting Point as Function of the Relative Roughnessε/D (When Re⟶∞)

The special case of the Colebrook equation when Re⟶∞ physically means that the flow friction factor λ in that case depends only on ε/D, for Re⟶∞, λ = f(ε/D), that is, the flow friction factor λ is not implicitly given [14]. In that way, the starting point can be calculated using the explicit equation which has only one variable, λ0 = f(ε/D) (Equation (2)). The results obtained in that way are accurate only for the case Re⟶∞ but for the smaller values of the Reynolds number Re which corresponds to the smooth turbulent flow, the error can goes up to 80% [13, 57]. Anyway, in that way, calculated value can be efficiently used as an initial starting guess for the iterative procedures for the whole domain of applicability of the Colebrook equation.(2)x0=1λ0=−2⋅log10ε/3.7⋅D︸roughpartofColebrook.The initial starting point obtained using the previous equation is referred as “traditional,” and 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 the accuracy ofλi+1 − λi ≤ 10−8, usually 6 steps are enough regarding the Newton–Raphson method (Figure 1).Figure 1
“Slow area” which requires 6 iterations to reach accuracy of 10−8 regarding the “traditional” option for the starting point calculated through Equation (2), for solving Colebrook’s equation in Newton’s procedure when calculation goes through the transmission factor x.


### 2.1.2. Starting Point Obtained Using Approximations to the Colebrook Equation

Every approximation to the Colebrook equation,λ ≈ f(Re, ε/D), can be used to put the initial starting point as close as possible near the final accurate solution [30]. For example, using one of the approximations with the error of up to 10% for calculation of the initial starting point x0=1/λ, the final accurate value of the flow friction factor λ is reached in the worst-case scenario after 3 iterations using the Colebrook equation and one of the procedures from Section 3.2. After 3 iterations, the whole practical domain of applicability of the Colebrook equation is covered with the difference between the two final iterations less than 10−8, λi+1 − λi ≤ 10−8 (Figure 2). In average, the method requires 2.7 iterations in average for all cases with the set precision (stopping criterion) very close to zero (about 10−8) when calculation goes through the transmission factor x. The results from Figure 1 are from the 65536 pairs of the Reynolds number Re and the relative roughness ε/D over the domain of applicability of the Colebrook equation domain (values of the Reynolds number Re between 4000 and 108 and the relative roughness ε/D between 0 and 0.05, dividing them into 256 points each).Figure 2
Area in which 2 iterations (left) and 3 iterations (right) are sufficient to calculate the flow friction factorλ with accuracy of 10−8 using Colebrook equation solved in Newton’s procedure when calculation goes through the transmission factor x and using approximation with error of up to 10%.


## 2.1.1. Starting Point as Function of the Relative Roughnessε/D (When Re⟶∞)

The special case of the Colebrook equation when Re⟶∞ physically means that the flow friction factor λ in that case depends only on ε/D, for Re⟶∞, λ = f(ε/D), that is, the flow friction factor λ is not implicitly given [14]. In that way, the starting point can be calculated using the explicit equation which has only one variable, λ0 = f(ε/D) (Equation (2)). The results obtained in that way are accurate only for the case Re⟶∞ but for the smaller values of the Reynolds number Re which corresponds to the smooth turbulent flow, the error can goes up to 80% [13, 57]. Anyway, in that way, calculated value can be efficiently used as an initial starting guess for the iterative procedures for the whole domain of applicability of the Colebrook equation.(2)x0=1λ0=−2⋅log10ε/3.7⋅D︸roughpartofColebrook.The initial starting point obtained using the previous equation is referred as “traditional,” and 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 the accuracy ofλi+1 − λi ≤ 10−8, usually 6 steps are enough regarding the Newton–Raphson method (Figure 1).Figure 1
“Slow area” which requires 6 iterations to reach accuracy of 10−8 regarding the “traditional” option for the starting point calculated through Equation (2), for solving Colebrook’s equation in Newton’s procedure when calculation goes through the transmission factor x.


## 2.1.2. Starting Point Obtained Using Approximations to the Colebrook Equation

Every approximation to the Colebrook equation,λ ≈ f(Re, ε/D), can be used to put the initial starting point as close as possible near the final accurate solution [30]. For example, using one of the approximations with the error of up to 10% for calculation of the initial starting point x0=1/λ, the final accurate value of the flow friction factor λ is reached in the worst-case scenario after 3 iterations using the Colebrook equation and one of the procedures from Section 3.2. After 3 iterations, the whole practical domain of applicability of the Colebrook equation is covered with the difference between the two final iterations less than 10−8, λi+1 − λi ≤ 10−8 (Figure 2). In average, the method requires 2.7 iterations in average for all cases with the set precision (stopping criterion) very close to zero (about 10−8) when calculation goes through the transmission factor x. The results from Figure 1 are from the 65536 pairs of the Reynolds number Re and the relative roughness ε/D over the domain of applicability of the Colebrook equation domain (values of the Reynolds number Re between 4000 and 108 and the relative roughness ε/D between 0 and 0.05, dividing them into 256 points each).Figure 2
Area in which 2 iterations (left) and 3 iterations (right) are sufficient to calculate the flow friction factorλ with accuracy of 10−8 using Colebrook equation solved in Newton’s procedure when calculation goes through the transmission factor x and using approximation with error of up to 10%.


## 2.2. Fixed Initial Starting Point

An idea from geometry to find “center of gravity” is used for the points for which the Newton–Raphson, the Halley, Schröder, and the 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 [58].


### 2.2.1. Fixed Initial Starting Point for the Newton–Raphson Method

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


### 2.2.2. Fixed Initial Starting Point for the Halley and Schröder Method

The starting point for calculation through the Halley and the Schröder method using Equation (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’s and Schröder’s method that requires 4 iterations using staring point Equation (2) has coordinates: (log10 (Re) = 5.3108⟶Re ≈ 204550; –log10 (ε/D) = 4.9431⟶ε/D ≈ 1.14·10−5)⟶λ0 = 0.015663210285978339, that is, x0 = 7.990256504. This value is the new optimal initial starting point in the case of the Halley and the Schröder method.Figure 4
“Slow area” which requires 4 iterations to reach accuracy of 10−8 regarding the “traditional” option for the starting point calculated through Equation (2), for solving Colebrook’s equation in Halley’s and Schröder’s procedure when calculation goes through the transmission factor x.With the new initial starting pointx0 = 7.990256504, three iterations are required at maximum to reach the required accuracy in case of the Halley and the Schröder method (Figure 5) as described in Section 3.2.Figure 5
Decreased maximal number of required iterations from 4 to 3 to reach accuracy of 10−8 for solving Colebrook’s equation in Halley’s and Schröder’s procedure when calculation goes through the transmission factor x where the initial starting point is with the fixed value: x0 = 7.990256504.


### 2.2.3. Fixed Initial Starting Point for the Three-Point Iterative Methods

The optimal normalized parameters for the fixed initial starting point for the three-point iterative methods explained in Section3.4. of this paper [26–28] are as follows: (log10 (Re) = 4.90060379974617⟶Re = 79543.33576; –log (ε/D)  =  5.33355157079189 ⟶ε/D  =  4.63926·10−6)⟶λ0 = 0.018904186734624⟶x0 = 7.273124147. The Džunić–Petković method is shown in Section 3.4. Additional recommendations about the initial starting point regarding the three-point iterative methods can be found in [59].


## 2.2.1. Fixed Initial Starting Point for the Newton–Raphson Method

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


## 2.2.2. Fixed Initial Starting Point for the Halley and Schröder Method

The starting point for calculation through the Halley and the Schröder method using Equation (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’s and Schröder’s method that requires 4 iterations using staring point Equation (2) has coordinates: (log10 (Re) = 5.3108⟶Re ≈ 204550; –log10 (ε/D) = 4.9431⟶ε/D ≈ 1.14·10−5)⟶λ0 = 0.015663210285978339, that is, x0 = 7.990256504. This value is the new optimal initial starting point in the case of the Halley and the Schröder method.Figure 4
“Slow area” which requires 4 iterations to reach accuracy of 10−8 regarding the “traditional” option for the starting point calculated through Equation (2), for solving Colebrook’s equation in Halley’s and Schröder’s procedure when calculation goes through the transmission factor x.With the new initial starting pointx0 = 7.990256504, three iterations are required at maximum to reach the required accuracy in case of the Halley and the Schröder method (Figure 5) as described in Section 3.2.Figure 5
Decreased maximal number of required iterations from 4 to 3 to reach accuracy of 10−8 for solving Colebrook’s equation in Halley’s and Schröder’s procedure when calculation goes through the transmission factor x where the initial starting point is with the fixed value: x0 = 7.990256504.


## 2.2.3. Fixed Initial Starting Point for the Three-Point Iterative Methods

The optimal normalized parameters for the fixed initial starting point for the three-point iterative methods explained in Section3.4. of this paper [26–28] are as follows: (log10 (Re) = 4.90060379974617⟶Re = 79543.33576; –log (ε/D)  =  5.33355157079189 ⟶ε/D  =  4.63926·10−6)⟶λ0 = 0.018904186734624⟶x0 = 7.273124147. The Džunić–Petković method is shown in Section 3.4. Additional recommendations about the initial starting point regarding the three-point iterative methods can be found in [59].


## 2.3. Starting Point for the Lambert W-Expressed Colebrook Equation

The friction factorλ in the Colebrook equation can be expressed in the explicit way through the Lambert W-function [3, 5, 8, 30, 42, 60, 61]. The Lambert W-function can further be evaluated using some types of the Householder iterative methods as shown in Section 3.5 of this paper.The Colebrook equation in a closed form through the Lambert W-function can be expressed in two ways, Equations (3) and (4). The first expression is [3, 30, 60, 62] as follows:(3)1λ=−2⋅log102⋅2.51⋅WyRe⋅ln10+ε3.7⋅D=−2⋅log1010−Wy/ln10+ε3.7⋅D,where y=Re⋅ln10/2⋅2.51≈Re/2.18.The argument of the Lambert W-function in this case depends only on the Reynolds number,λ0 = Re/2.18. Knowing that the practical range of the Reynolds number goes from 4000 to 108, the argument of the Lambert W-function is goes from about 1835 to 45871560, where W(1835) = 5.763291081 and W(45871560) = 14.93748223. The Halley procedure is fast and any initial starting point can be chosen between 5 and 15, but the Newton–Raphson method is very slow and we found that, for the best results, the initial starting point 15 has to be chosen. Note that for W(45871560) = 14.93748223, the Newton–Raphson procedure does not work in Excel for the values of initial starting point lower than 8.814.Due to transformations of coefficients, Equation (3) can introduce the relative error up to 2% and should be considered as explicit approximation to the Colebrook equation rather to its equivalent [63].The second expression is as follows [42, 54]:(4)λ=2ln10⋅Weα−Re⋅ε/D3.7⋅2.51−2,where α=Re⋅ε/D⋅ln10/2⋅2.51⋅3.7−ln2⋅3.7/Re⋅ln10.Argument of the Lambert W-function in this case iseα which for the certain combinations of the Reynolds number Re and the relative roughness ε/D from the practical domain of the Colebrook equation is too big to be calculated in registers of computers [6, 41, 54]. This can be overwhelmed with the Wright Ω-function, ωβ=Weα [40, 43, 64–66].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 [41, 54, 67, 68].


## 3. Iterative Methods Adopted for the Colebrook Equation

The Householder method [22] is a numerical algorithm for solving the nonlinear equation such as Colebrook’s. During the Householder procedure, in successive calculation, that is, in iterative cycles, the original assumed value of the unknown quantity (the initial starting point [50]) 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 [28].The following types of the method are used in this paper: the first-order Householder method (Newton–Raphson [18, 19]), the 2nd order (Halley [69, 70] and Schröder [25]), and the 3rd order, as well as the three-point methods [28]. All these methods require the calculation of the derivatives which is usually underlined as the most important shortcoming of the Householder and the three-point methods compared with the simple fixed-point procedure in respect of the Colebrook equation. The Newton–Raphson and the three-point methods (in most cases) require only the first derivative, the Halley and the Schröder method requires the first and the second derivative, while the 3rd 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 were 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 variant of the Newton–Raphson method [71].All shown approaches with the Householder methods in our case usually require only 2 to 4 iterations to reach the final accurate solution [72]. This number can be slightly higher depending on the chosen method where the secant method requires by default 1-2 iterations more. Also some simple transformations of the Colebrook equation, such as introduction of the transmission factor in form of the shift x=1/λ, can reduce the number of required iterations. Knowing that the right form of equation is essential for all types of the Householder methods, here are examined the two at first look very similar options: (1) direct calculation of λ in Section 3.1 and (2) indirect calculation of λ through transmission factor x=1/λ in Section 3.2.Finally, the Colebrook equation can be rewritten in an explicit form through the Lambert W-function [41, 42, 54, 73] and the Lambert W-function is solved in Section 3.5 using the Newton–Raphson and the Halley procedure.


### 3.1. Direct Calculation ofλ with Derivative Calculated in Analytical Way

The proposed technique requires the Colebrook equation in the formf(λ, Re, ε/D) = 0, Equation (5), where λ is treated as variable, the first derivative f′(λ, Re, ε/D), Equation (6) of the Colebrook equation with respect to λ 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 formf(λ, Re, ε/D) = 0:(5)fλ=1/λ+2⋅log102.51/Re⋅λ+ε/3.7⋅D=0︸Colebrook:fλ=0.The first derivativef′ with respect to λ in exact analytical way:(6)f′λ=ddλfλ=−1/2⋅1/λ3⋅1+2⋅2.51/ln10⋅Re⋅2.51/Re⋅λ+ε/3.7⋅D︸1st  derivativef′λ−analytical.Initial valueλ0 is selected as explained in Section 2 of this paper in order to calculate the residue f/f′ and start the iterative procedure:(7)λ1=λ0−fλ0/f′λ0︸Newton–Raphson=λ0−1/λ0+2⋅log102.51/Re⋅λ0+ε/3.7⋅D−1/2⋅1/λ03⋅1+2.18/Re⋅2.51/Re⋅λ0+ε/3.7⋅D.The procedureλi+1 = λi − f(λi)/f′(λi) needs to be followed until the residue f(λi)/f′(λi) ≈ 0.The explained Newton–Raphson procedure is shown in Tables1 and 2 for two numerical examples: (1) (Re = 5·106, ε/D = 2.5·10−5)⟶λ = 0.010279663295529 and (2) (Re = 3·104, ε/D = 9·10−3)⟶λ = 0.038630738574792. As explained in Section 2 of this paper, the initial starting point λ0 in Table 1 depends on the input parameters, while in Table 2 it is with the fixed value.Table 1
Newton–Raphson procedure. Option 1: starting pointλ0 depends on input parameters: Equation (2), calculation of λ: Equation (7), and analytical derivative f′(λ): Equation (6).
Re = 5·106, ε/D = 2.5·10−5
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.009352225155363
Iteration 1
0.495092014
−573.0134258
0.010216239839661
Iteration 2
0.031705666
−502.2190127
0.010279370993451
Iteration 3
0.000145453
−497.622807
0.010279663289327
Iteration 4
0.000000003
−497.6016902
λ = 0.010279663295529
Control step
0.000000000
−497.6016898
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.036588313752304
Iteration 1
0.143632267
−73.25157738
0.038549121591193
Iteration 2
0.005520057
−67.74092562
0.038630609361351
Iteration 3
0.000008725
−67.52696208
0.038630738574469
Iteration 4
0.000000000
−67.5266237
λ = 0.038630738574792
Control step
0.000000000
−67.5266237
0.038630738574792Table 2
Newton–Raphson procedure. Option 2: fixed initial starting pointλ0 = 0.024069128765100981 from Section 2.2.1, calculation of λ: Equation (7), and analytical derivative f′(λ): Equation (6).
Re = 5·106, ε/D = 2.5·10−5
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.024069128765101
Iteration 1
−3.554956084
−139.7424853
−0.001370207567104
Iteration 2
17.630891548
−10069.59089
0.000380696888310
Iteration 3
42.275315189
−68216.8306
0.001000416608714
Iteration 4
22.325487096
−16105.99979
0.002386576262278
Iteration 5
10.932300910
−4398.30144
0.004872149626988
Iteration 6
4.615550920
−1516.202309
0.007916302041016
Iteration 7
1.426053458
−734.846953
0.009856914916156
Iteration 8
0.217044469
−529.7853757
0.010266598684182
Iteration 9
0.006507144
−498.5470019
0.010279650902858
Iteration 10
0.000006167
−497.602585
0.010279663295518
Iteration 11
0.000000000
−497.6016898
λ = 0.010279663295529
Control step
0.000000000
−497.6016898
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.024069128765101
Iteration 1
1.391712394
−137.1740994
0.034214720386916
Iteration 2
0.326434508
−80.9945153
0.038245048943635
Iteration 3
0.026240732
−68.54940037
0.038627849256271
Iteration 4
0.000195117
−67.53419088
0.038630738412914
Iteration 5
0.000000011
−67.52662412
λ = 0.038630738574792
Control step
0.000000000
−67.5266237
0.038630738574792Here shown direct calculation of the unknown flow friction factorλ is sensitive on the chosen initial starting point λ0 [50]. The fixed initial point λ0 chosen as in Section 2.2 in some cases requires the increased number of iterations to reach the final solution although the procedure still maintains very good convergent properties [71, 72]. To reduce the number of the required iterations, use of some of the explicit approximations to the Colebrook equation is advised in order to bring the initial starting point λ0 as close as possible near the final calculated value. Therefore, the approach with the fixed starting point as explained in Section 2.2 of this paper, in this case, cannot be advised in comparison with the approach with the starting point obtained using approximations as explained in Section 2.1.Comparing the same approach but with the different starting points (Tables1 and 2), we can conclude that the one single calculated negative value for flow friction factor λ can increase the 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.


### 3.2. Indirect Calculation ofλ through the Transmission Factor x=1/λ

The appropriate form of the function is essential to reduce the number of required iteration to reach the final solution. In order to accelerate the procedure, an appropriate shiftx=1/λ is used to provide some kind of linearization of the problem.The Newton–Raphson procedure with these changes has similar steps as already shown:Shift in form of the transmission factorx=1/λ should be introduced in order to transform the Colebrook equation in form f(x, Re, ε/D) = 0:(8)fx=x+2⋅log102.51⋅x/Re+ε/3.7⋅D=0︸Colebrook:fx=0.The first derivative of Equation (8) in respect of 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.


#### 3.2.1. Indirect Calculation ofλ through the Transmission Factor x=1/λ with the Derivative Calculated Analytically

The first derivativef′ with respect to x can be obtained analytically, Equation (9) (also Equation (11) gives the same results):(9)f′x=ddxfx=1+2⋅2.51/Re⋅ln10/ε/3.7⋅D+2.51/Re⋅x︸1st  derivativef′x−analytical.Initial value of the flow friction factorλ0 should be chosen and the residue f/f′ is calculated in order to start the Newton–Raphson procedure:(10)x1=x0−fx0/f′x0︸Newton−Raphson=x0−x0+2⋅log102.51⋅x0/Re+ε/3.7⋅D1+2⋅2.51/Re⋅ln10/ε/3.7⋅D+2.51/Re⋅x0.The procedurexi+1 = xi − f(xi)/f′(xi) should be followed until the residue f(xi)/f′(xi) ≈ 0. Then, the final solution is λn=xn−2, 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 x0 = 6.445695939 for Re = 5·106, ε/D = 2.5·10−5).Table 3
Newton–Raphson procedure. Option 3: fixed initial starting pointx0 = 6.445695939→λ0 = 0.024069128765101 from Section 2.2.1, indirect calculation of λ through the transmission factor x: Equation (10), and analytical derivative f′(x): Equation (9).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (9)
x
0 = 6.445695939
λ
0 = 0.024069128768719
Iteration 1
−3.554956085
1.043635910
9.852014225862620
0.010302673560706
Iteration 2
−0.011430857
1.037259804
9.863034470914730
0.010279663490514
Iteration 3
−0.000000097
1.037242198
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (9)
x
0 = 6.445695939
λ
0 = 0.024069128768719
Iteration 1
1.391712393
1.024454486
5.087204750239650
0.038640395682209
Iteration 2
−0.000651990
1.025427001
5.087840572945700
0.038630738577020
Iteration 3
0.000000000
1.025426528
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.046830475
5.087840573092420
0.038630738574792


#### 3.2.2. Indirect Calculation ofλ through the Transmission Factor x=1/λ with the Symbolic Derivative

The exact analytical expression of the first derivativef′ with respect to x can be obtained in MATLAB, Equation (11); results are the same as using Equation (9):(11)f′x=ddxfx=5.02/Re⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re+1=9287⋅ln10⋅x+1000⋅ln10⋅ε/D⋅Re+18574/ln10⋅9287⋅x+1000⋅ε/D⋅Re︸1st  derivativef′x−MATLAB.Initial value of the flow friction factorλ0 should be chosen and the residue f/f′ is calculated in order to start the Newton–Raphson procedure:(12)x1=x0−fx0/f′x0︸Newton−Raphson=x0−x0+2⋅log102.51⋅x0/Re+ε/3.7⋅D21384.11⋅x0+2302.58⋅ε/D⋅Re+18574/21384.11⋅x0+2302.58⋅ε/D⋅Re.The procedurexi+1 = xi − f(xi)/f′(xi) should be followed until the residue f(xi)/f′(xi) ≈ 0. Then, the final solution is λn=xn−2, where n = i + 1 is the final iteration.The iterative procedure can be accelerated using Halley’s formula instead of the Newton–Raphson:(13)x1=x0−fx0/f′x0/1−fx0·f″x0/2⋅f′x02︸Halley=x0−2⋅fx0⋅f′x0/2⋅f′x02−fx0⋅f″x0︸Halley.In general,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is the final iteration in which xn ≈ xn+1.The second derivativef″(x) in respect of x is required:(14)f″x=ddxf′x=−12.6/Re2⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re2=−172496738/ln10⋅9287⋅x+1000⋅ε/D⋅Re2︸2ndderivativef″x−MATLAB.The Newton–Raphson method belongs to the 1st order and the Halley to the 2nd order of Householder’s method, while the 3rd order can be expressed using the following equation:(15)x1=x0−6⋅fx0⋅f′x02−3⋅fx02⋅f″x0/6⋅f′x03−6⋅fx0⋅f′x0⋅f″x0+fx02⋅f‴x0︸3rdorderHouseholder.Again,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is the final iteration in which xn ≈ xn+1.The required 3rd derivative f‴(x) can be expressed using the following equation:(16)f‴x=ddxf″x=63.253/Re3⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re3=32039544116.2/ln10⋅9287⋅x+1000⋅ε/D⋅Re3︸3rd  derivativef‴x−MATLAB.Also here one has to be underlined that the Halley method [74] is not the unique Householder method of the 2nd order [22]. For example, the Schröder method [25] also belongs to the group:(17)x1=x0−fx0/f′x0−f″x0⋅fx02/2⋅f′x03︸Schröder.Further,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is final iteration in which xn ≈ xn+1.Using the presented Householder procedures, the 1st order: the Newton–Raphson, the 2nd order: Halley, and the 3rd order, the unknown flow friction factor λ should be calculated for the two given pairs of the Reynolds number Re and the relative roughness ε/D: (1) (Re = 5·106, ε/D = 2.5·10−5)⟶λ = 0.010279663295529 and (2) (Re = 3·104, ε/D = 9·10−3)⟶λ = 0.038630738574792. The calculation presented in Tables 4–7 is through the transmission factor x, using the symbolic derivative f′(x), but with the different initial starting point λ0.Table 4
Newton–Raphson procedure. Option 4: starting pointx0 depends on input parameters: Equation (2), indirect calculation of λ through the transmission factor x: Equation (12), and the symbolic derivative f′(x): Equation (11).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.495092014
1.036495031
9.862863625818000
0.010280019623455
Iteration 2
−0.000177305
1.037242471
9.863034564433310
0.010279663295576
Iteration 3
0.000000000
1.037242198
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
x
0 = 5.227918429
λ
0 = 0.036588313752304
Iteration 1
0.143632267
1.025322691
5.087833489750430
0.038630846139210
Iteration 2
−0.000007263
1.025426533
5.087840573092400
0.038630738574793
Iteration 3
0.000000000
1.025426528
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
5.087840573092420
0.038630738574792Table 5
Halley procedure. Option 5: fixed initial starting pointx0 = 7.990256504→λ0 = 0.015663210285978339 from Section 2.2.2, indirect calculation of λ through the transmission factor x: Equation (13), and the symbolic derivatives f′(x) and f″(x): Equations (11) and (14).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
−1.945484250
1.040493788
−0.001887828
9.863203600915390
0.010279310950983
Iteration 2
0.000175332
1.037241928
−0.001596798
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
2.973246188
1.023435376
−0.000632309
5.087698791122220
0.038632891696967
Iteration 2
−0.000145387
1.025426633
−0.000744326
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
5.087840573092420
0.038630738574792Table 6
Third-order Householder’s procedure. Option 6: starting pointx0 depends on input parameters: Equation (2), indirect calculation of λ through the transmission factor x: Equation (15), and the symbolic derivatives f′(x), f″(x), and f‴(x): Equations (11), (14), and (16).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
f‴(x), Equation (16)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.495092014
1.036495031
−0.001533392
0.000128855
9.863034531578420
0.010279663364062
Iteration 2
−0.000000034
1.037242198
−0.001596821
0.000136933
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
0.000136933
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
f‴(x), Equation (16)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.143632267
1.025322691
−0.000738253
0.000043046
5.087840573035260
0.038630738575660
Iteration 2
0.000000000
1.025426528
−0.000744320
0.000043578
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
0.000043578
5.087840573092420
0.038630738574792Table 7
Schröder procedure. Option 7: fixed initial starting pointx0 = 7.990256504→λ0 = 0.015663210285978339 from Section 2.2.2, indirect calculation of λ through the transmission factor x: Equation (17), and the symbolic derivatives f′(x) and f″(x): Equations (11) and (14).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
−1.945484250
1.040493788
−0.001887828
9.863198212166060
0.010279322183170
Iteration 2
0.000169742
1.037241937
−0.001596799
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
2.973246188
1.023435376
−0.000632309
5.087701128882780
0.038632856193927
Iteration 2
−0.000142990
1.025426632
−0.000744326
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
5.087840573092420
0.038630738574792


### 3.3. Secant Method

The secant method is similar to the Newton–Raphson; it requires two starting pointsλ0 and λ−1 but does not require calculation of the derivatives [71]. The approach with the direct calculation of λ with the two required starting points λ0 and λ−1 is formulated:(18)λ1=λ0−fλ0fλ−1−fλ0/λ−1−λ0.Counteri starts from i − 1 and goes to n + 1 in which λn = λn+1.The approach through the transmission factor x is as follows:(19)x1=x0−fx0fx−1−fx0/x−1−x0.As already described, counteri also starts from i − 1 and goes to n + 1 in which xn = xn+1.The flow friction factorλ is calculated in Tables 8 and 9 for two pairs of the Reynolds number and the relative roughness (1) Re = 5·106, ε/D = 2.5·10−5, and (2) Re = 3·104, ε/D = 9·10−3, using the secant procedure with direct calculation of λ and indirect through the transmission factor x.Table 8
Secant procedure. Option 8: two initial starting pointsλ0 and λ−1 required: starting point λ−1 is with fixed value λ−1 = 0.024069128765101 (i.e., x−1 = 6.445695939) as in Section 2.2.1, while starting point λ0 depends on input parameters: Equation (2), and direct calculation of λ: Equation (18).
Re = 5·106, ε/D = 2.5·10−5
f(λi−1), Equation (5)
f(λi), Equation (5)
f
λ
i
−
1
−
f
λ
i
/
λ
i
−
1
−
λ
i
λ
−1 = 0.024069128765101
λ
0 = 0.009352225155363
Iteration 1
0.495092014
−3.554956084
−275.1970255
0.011151270814558
Iteration 2
−0.408071981
0.495092014
−502.0239429
0.010338417191085
Iteration 3
−0.029111936
−0.408071981
−466.2094551
0.010275973292109
Iteration 4
0.001836644
−0.029111936
−495.6221591
0.010279679026163
Iteration 5
−0.000007828
0.001836644
−497.734448
0.010279663299743
Iteration 6
−0.000000002
−0.000007828
−497.6011214
λ = 0.010279663295529
Control step
0.000000000
−0.000000002
−497.6012179
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(λi−1), Equation (5)
f(λi), Equation (5)
f
λ
i
−
1
−
f
λ
i
/
λ
i
−
1
−
λ
i
λ
−1 = 0.024069128765101
λ
0 = 0.036588313752304
Iteration 1
1.391712394
0.143632267
−99.69340079
0.038029053721052
Iteration 2
0.143632267
0.041110009
−71.15944585
0.038606770549177
Iteration 3
0.041110009
0.001619232
−68.35663251
0.038630458556837
Iteration 4
0.001619232
0.000018909
−67.5583902
0.038630738444645
Iteration 5
0.000018909
0.000000009
−67.52699052
λ = 0.038630738574792
Control step
0.000000009
0.000000000
−67.52662212
0.038630738574792Table 9
Secant procedure. Option 9: two initial starting pointsx0 and x−1 required: starting point λ−1 is with fixed value x−1 = 6.445695939 (i.e., λ−1 = 0.024069128765101) as in Section 2.2.1, while starting point λ0 depends on input parameters: Equation (2), and indirect calculation of λ through the transmission factor x: Equation (19).
Re = 5·106, ε/D = 2.5·10−5
f(xi−1), Equation (8)
f(xi), Equation (8)
f
x
i
−
1
−
f
x
i
/
x
i
−
1
−
x
i
x
−1 = 6.445695939
λ
−1 = 0.024069128765101
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
−3.554956084
0.495092014
1.039853012
9.864406125318800
0.010276804896656
Iteration 2
0.495092014
0.001422639
1.03686501
9.863034066961850
0.010279664332547
Iteration 3
0.001422639
−0.000000516
1.037241104
9.863034564456330
0.010279663295528
Iteration 4
−0.000000516
0.000000000
1.0372422
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
0.000000000
1.037162162
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(xi−1), Equation (8)
f(xi), Equation (8)
f
x
i
−
1
−
f
x
i
/
x
i
−
1
−
x
i
x
−1 = 6.445695939
λ
−1 = 0.024069128765101
x
0 = 5.227918429
λ
0 = 0.036588313752304
Iteration 1
1.391712394
0.143632267
1.024883541
5.087773465040530
0.038631757665255
Iteration 2
0.143632267
−0.000068814
1.025374564
5.087840576494990
0.038630738523123
Iteration 3
−0.000068814
0.000000003
1.025426553
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000003
0.000000000
1.025426591
5.087840573092420
0.038630738574792Calculation 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·106, ε/D = 2.5·10−5, and the required number of iterations is 4 in direct calculation and 3 in indirect for Re = 3·104, ε/D = 9·10−3.


### 3.4. Three-Point Methods

Three-point iterative methods require in every iteration an evaluation of the functionf at three points: x0, y0, and z0. However, these methods converge very fast to the accurate solution. The mathematical background of the three-point iterative methods is given by Sharma and Arora [28]. Here, we will apply the Džunić–Petković–Petković three-point iterative method to the Colebrook equation [26, 28]. It requires only one iteration (up to 2 in the worst cases) to reach the final accurate solution and we will show all steps to calculate the friction factor λ, Equation (20), where the numerical values are for Re = 5·106, ε/D = 2.5·10−5. The initial starting point is x0 = 7.273124147 as described in Section 2.2.3.(20)x0=7.273124147,fx0=x0+2⋅log102.51⋅x0Re+ε3.7⋅D=−2.692152546,f′x0=9287⋅ln10⋅x0+1000⋅ln10⋅ε/D⋅Re+18574ln10⋅9287⋅x0+1000⋅ε/D⋅Re=5.02Re⋅ln10⋅10/37⋅ε/D+2.51⋅x0/Re+1=1.041894438,y0=x0−fx0f′x0=7.273124147−−2.6921525461.041894=9.857025593360860,fy0=y0+2⋅log102.51⋅y0Re+ε3.7⋅D=−0.006232787,z0=y0−fx0fx0−2·fy0·fy0f′x0=9.863035589,fz0=z0+2⋅log102.51⋅z0Re+ε3.7⋅D=−0.006232787,x1=z0−fz0f′x0·1−2·fy0/fx0−fy0/fx02·1−fz0/fy0·1−2·fz0/fx0=9.863034564,x1=9.863034564⟶λ1=0.010279663295529.Three-point methods are also shown in detail in Praks and Brkić [75].


### 3.5. Expressed through the Lambert W-Function

The mathematical background of the Lambert W-function is presented in the literature [44–46]. Moreover, an application of the Lambert W-function in hydraulic problems is given by [41, 54, 73].The Lambert W-functionW(y) [46] is the solution of y=z⋅ez, which needs to be in the appropriate form:(21)fz=z⋅ez−y=0.The first derivativef′(z) is as follows:(22)f′z=ez⋅z+1.Choose initial valuez0, calculate the residue f/f′, and start the procedure:(23)z1=z0−fz0f′z0=z0−z0⋅ez0−y/ez0⋅z0+1︸Newton−Raphson=z0−z0⋅ez0−Re/2.18ez0⋅z0+1.Then follow the procedurezi+1 = zi − f(zi)/f′(zi) until the residue f(zi)/f′(zi)≈0, where n = i + 1 in final iteration.Halley’s procedure:(24)z1=z0−fz0f′z0−fz0⋅f″z0/2⋅f′z0=z0−z0⋅ez0−y/ez0⋅z0+1−z0⋅ez0−y⋅z0+2/2⋅z0+1︸Halley,where the second derivative is(25)f′z=ez⋅z+2.The Schröder expression is as follows:(26)z1=z0−fz0f′z0−f″z0⋅fz022⋅f′z03=z0−z⋅ez−y/ez⋅z+1−ez⋅z+2⋅z⋅ez−y2/2⋅ez⋅z+13︸Schröder.Further, in all cases, the Newton–Raphson, the Halley, and the Schröder;z1 = zi+1 and z0 = zi; i = 0 to n, where n + 1 is final iteration in which zn ≈ zn+1.The argument of the Lambert W-functiony, in our case defined by Equation (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 the Halley method for Re = 5·106 and Re = 3·104 where initial starting point is set as z0 = 15 as recommended in Section 2.3 of this paper (for z0 < 8.814, the Newton–Raphson procedure cannot start).Table 10
Calculation of W(y) wherey = Re/2.18 using the Newton–Raphson, the Halley, and the Schröder iterative methods.
Re = 5·106, y = Re/2.18 = 2293411.45
Re = 3·104, y = Re/2.18 = 13760.47
Newton–Raphson
Halley
Schröder
Newton–Raphson
Halley
Schröder
Iteration 0
z
0 = 15
z
0 = 15
z
0 = 15
z
0 = 15
z
0 = 15
z
0 = 15
Iteration 1
14.10634749
13.29860556
13.68208338
14.06276308
13.1333396
13.59610616
Iteration 2
13.28604947
12.2757343
12.62556802
13.12986539
11.29171838
12.20340124
Iteration 3
12.62863905
12.14855784
12.17738057
12.20257063
9.520829163
10.83006437
Iteration 4
12.25343232
z = 12.14835704
12.14836628
11.28354302
8.068323472
9.505729616
Iteration 5
12.15407754
12.14835704
z = 12.14835704
10.37904335
7.530266826
8.341483562
Iteration 6
12.14837461
—
12.14835704
9.504505014
7.512930233
7.637280252
Iteration 7
z = 12.14835704
—
—
8.697314341
z = 7.512929679
7.513654122
Iteration 8
12.14835704
—
—
8.037456295
7.512929679
z = 7.512929679
Iteration 9
—
—
—
7.640105762
—
7.512929679
Iteration 10
—
—
—
7.52154464
—
—
Iteration 11
—
—
—
7.512971011
—
—
Iteration 12
—
—
—
7.51292968
—
—
Iteration 13
—
—
—
z = 7.512929679
—
—
Iteration 14
—
—
—
7.512929679
—
—


## 3.1. Direct Calculation ofλ with Derivative Calculated in Analytical Way

The proposed technique requires the Colebrook equation in the formf(λ, Re, ε/D) = 0, Equation (5), where λ is treated as variable, the first derivative f′(λ, Re, ε/D), Equation (6) of the Colebrook equation with respect to λ 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 formf(λ, Re, ε/D) = 0:(5)fλ=1/λ+2⋅log102.51/Re⋅λ+ε/3.7⋅D=0︸Colebrook:fλ=0.The first derivativef′ with respect to λ in exact analytical way:(6)f′λ=ddλfλ=−1/2⋅1/λ3⋅1+2⋅2.51/ln10⋅Re⋅2.51/Re⋅λ+ε/3.7⋅D︸1st  derivativef′λ−analytical.Initial valueλ0 is selected as explained in Section 2 of this paper in order to calculate the residue f/f′ and start the iterative procedure:(7)λ1=λ0−fλ0/f′λ0︸Newton–Raphson=λ0−1/λ0+2⋅log102.51/Re⋅λ0+ε/3.7⋅D−1/2⋅1/λ03⋅1+2.18/Re⋅2.51/Re⋅λ0+ε/3.7⋅D.The procedureλi+1 = λi − f(λi)/f′(λi) needs to be followed until the residue f(λi)/f′(λi) ≈ 0.The explained Newton–Raphson procedure is shown in Tables1 and 2 for two numerical examples: (1) (Re = 5·106, ε/D = 2.5·10−5)⟶λ = 0.010279663295529 and (2) (Re = 3·104, ε/D = 9·10−3)⟶λ = 0.038630738574792. As explained in Section 2 of this paper, the initial starting point λ0 in Table 1 depends on the input parameters, while in Table 2 it is with the fixed value.Table 1
Newton–Raphson procedure. Option 1: starting pointλ0 depends on input parameters: Equation (2), calculation of λ: Equation (7), and analytical derivative f′(λ): Equation (6).
Re = 5·106, ε/D = 2.5·10−5
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.009352225155363
Iteration 1
0.495092014
−573.0134258
0.010216239839661
Iteration 2
0.031705666
−502.2190127
0.010279370993451
Iteration 3
0.000145453
−497.622807
0.010279663289327
Iteration 4
0.000000003
−497.6016902
λ = 0.010279663295529
Control step
0.000000000
−497.6016898
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.036588313752304
Iteration 1
0.143632267
−73.25157738
0.038549121591193
Iteration 2
0.005520057
−67.74092562
0.038630609361351
Iteration 3
0.000008725
−67.52696208
0.038630738574469
Iteration 4
0.000000000
−67.5266237
λ = 0.038630738574792
Control step
0.000000000
−67.5266237
0.038630738574792Table 2
Newton–Raphson procedure. Option 2: fixed initial starting pointλ0 = 0.024069128765100981 from Section 2.2.1, calculation of λ: Equation (7), and analytical derivative f′(λ): Equation (6).
Re = 5·106, ε/D = 2.5·10−5
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.024069128765101
Iteration 1
−3.554956084
−139.7424853
−0.001370207567104
Iteration 2
17.630891548
−10069.59089
0.000380696888310
Iteration 3
42.275315189
−68216.8306
0.001000416608714
Iteration 4
22.325487096
−16105.99979
0.002386576262278
Iteration 5
10.932300910
−4398.30144
0.004872149626988
Iteration 6
4.615550920
−1516.202309
0.007916302041016
Iteration 7
1.426053458
−734.846953
0.009856914916156
Iteration 8
0.217044469
−529.7853757
0.010266598684182
Iteration 9
0.006507144
−498.5470019
0.010279650902858
Iteration 10
0.000006167
−497.602585
0.010279663295518
Iteration 11
0.000000000
−497.6016898
λ = 0.010279663295529
Control step
0.000000000
−497.6016898
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(λ), Equation (5)
f′(λ), Equation (6)
λ
0 = 0.024069128765101
Iteration 1
1.391712394
−137.1740994
0.034214720386916
Iteration 2
0.326434508
−80.9945153
0.038245048943635
Iteration 3
0.026240732
−68.54940037
0.038627849256271
Iteration 4
0.000195117
−67.53419088
0.038630738412914
Iteration 5
0.000000011
−67.52662412
λ = 0.038630738574792
Control step
0.000000000
−67.5266237
0.038630738574792Here shown direct calculation of the unknown flow friction factorλ is sensitive on the chosen initial starting point λ0 [50]. The fixed initial point λ0 chosen as in Section 2.2 in some cases requires the increased number of iterations to reach the final solution although the procedure still maintains very good convergent properties [71, 72]. To reduce the number of the required iterations, use of some of the explicit approximations to the Colebrook equation is advised in order to bring the initial starting point λ0 as close as possible near the final calculated value. Therefore, the approach with the fixed starting point as explained in Section 2.2 of this paper, in this case, cannot be advised in comparison with the approach with the starting point obtained using approximations as explained in Section 2.1.Comparing the same approach but with the different starting points (Tables1 and 2), we can conclude that the one single calculated negative value for flow friction factor λ can increase the 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.


## 3.2. Indirect Calculation ofλ through the Transmission Factor x=1/λ

The appropriate form of the function is essential to reduce the number of required iteration to reach the final solution. In order to accelerate the procedure, an appropriate shiftx=1/λ is used to provide some kind of linearization of the problem.The Newton–Raphson procedure with these changes has similar steps as already shown:Shift in form of the transmission factorx=1/λ should be introduced in order to transform the Colebrook equation in form f(x, Re, ε/D) = 0:(8)fx=x+2⋅log102.51⋅x/Re+ε/3.7⋅D=0︸Colebrook:fx=0.The first derivative of Equation (8) in respect of 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.


### 3.2.1. Indirect Calculation ofλ through the Transmission Factor x=1/λ with the Derivative Calculated Analytically

The first derivativef′ with respect to x can be obtained analytically, Equation (9) (also Equation (11) gives the same results):(9)f′x=ddxfx=1+2⋅2.51/Re⋅ln10/ε/3.7⋅D+2.51/Re⋅x︸1st  derivativef′x−analytical.Initial value of the flow friction factorλ0 should be chosen and the residue f/f′ is calculated in order to start the Newton–Raphson procedure:(10)x1=x0−fx0/f′x0︸Newton−Raphson=x0−x0+2⋅log102.51⋅x0/Re+ε/3.7⋅D1+2⋅2.51/Re⋅ln10/ε/3.7⋅D+2.51/Re⋅x0.The procedurexi+1 = xi − f(xi)/f′(xi) should be followed until the residue f(xi)/f′(xi) ≈ 0. Then, the final solution is λn=xn−2, 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 x0 = 6.445695939 for Re = 5·106, ε/D = 2.5·10−5).Table 3
Newton–Raphson procedure. Option 3: fixed initial starting pointx0 = 6.445695939→λ0 = 0.024069128765101 from Section 2.2.1, indirect calculation of λ through the transmission factor x: Equation (10), and analytical derivative f′(x): Equation (9).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (9)
x
0 = 6.445695939
λ
0 = 0.024069128768719
Iteration 1
−3.554956085
1.043635910
9.852014225862620
0.010302673560706
Iteration 2
−0.011430857
1.037259804
9.863034470914730
0.010279663490514
Iteration 3
−0.000000097
1.037242198
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (9)
x
0 = 6.445695939
λ
0 = 0.024069128768719
Iteration 1
1.391712393
1.024454486
5.087204750239650
0.038640395682209
Iteration 2
−0.000651990
1.025427001
5.087840572945700
0.038630738577020
Iteration 3
0.000000000
1.025426528
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.046830475
5.087840573092420
0.038630738574792


### 3.2.2. Indirect Calculation ofλ through the Transmission Factor x=1/λ with the Symbolic Derivative

The exact analytical expression of the first derivativef′ with respect to x can be obtained in MATLAB, Equation (11); results are the same as using Equation (9):(11)f′x=ddxfx=5.02/Re⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re+1=9287⋅ln10⋅x+1000⋅ln10⋅ε/D⋅Re+18574/ln10⋅9287⋅x+1000⋅ε/D⋅Re︸1st  derivativef′x−MATLAB.Initial value of the flow friction factorλ0 should be chosen and the residue f/f′ is calculated in order to start the Newton–Raphson procedure:(12)x1=x0−fx0/f′x0︸Newton−Raphson=x0−x0+2⋅log102.51⋅x0/Re+ε/3.7⋅D21384.11⋅x0+2302.58⋅ε/D⋅Re+18574/21384.11⋅x0+2302.58⋅ε/D⋅Re.The procedurexi+1 = xi − f(xi)/f′(xi) should be followed until the residue f(xi)/f′(xi) ≈ 0. Then, the final solution is λn=xn−2, where n = i + 1 is the final iteration.The iterative procedure can be accelerated using Halley’s formula instead of the Newton–Raphson:(13)x1=x0−fx0/f′x0/1−fx0·f″x0/2⋅f′x02︸Halley=x0−2⋅fx0⋅f′x0/2⋅f′x02−fx0⋅f″x0︸Halley.In general,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is the final iteration in which xn ≈ xn+1.The second derivativef″(x) in respect of x is required:(14)f″x=ddxf′x=−12.6/Re2⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re2=−172496738/ln10⋅9287⋅x+1000⋅ε/D⋅Re2︸2ndderivativef″x−MATLAB.The Newton–Raphson method belongs to the 1st order and the Halley to the 2nd order of Householder’s method, while the 3rd order can be expressed using the following equation:(15)x1=x0−6⋅fx0⋅f′x02−3⋅fx02⋅f″x0/6⋅f′x03−6⋅fx0⋅f′x0⋅f″x0+fx02⋅f‴x0︸3rdorderHouseholder.Again,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is the final iteration in which xn ≈ xn+1.The required 3rd derivative f‴(x) can be expressed using the following equation:(16)f‴x=ddxf″x=63.253/Re3⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re3=32039544116.2/ln10⋅9287⋅x+1000⋅ε/D⋅Re3︸3rd  derivativef‴x−MATLAB.Also here one has to be underlined that the Halley method [74] is not the unique Householder method of the 2nd order [22]. For example, the Schröder method [25] also belongs to the group:(17)x1=x0−fx0/f′x0−f″x0⋅fx02/2⋅f′x03︸Schröder.Further,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is final iteration in which xn ≈ xn+1.Using the presented Householder procedures, the 1st order: the Newton–Raphson, the 2nd order: Halley, and the 3rd order, the unknown flow friction factor λ should be calculated for the two given pairs of the Reynolds number Re and the relative roughness ε/D: (1) (Re = 5·106, ε/D = 2.5·10−5)⟶λ = 0.010279663295529 and (2) (Re = 3·104, ε/D = 9·10−3)⟶λ = 0.038630738574792. The calculation presented in Tables 4–7 is through the transmission factor x, using the symbolic derivative f′(x), but with the different initial starting point λ0.Table 4
Newton–Raphson procedure. Option 4: starting pointx0 depends on input parameters: Equation (2), indirect calculation of λ through the transmission factor x: Equation (12), and the symbolic derivative f′(x): Equation (11).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.495092014
1.036495031
9.862863625818000
0.010280019623455
Iteration 2
−0.000177305
1.037242471
9.863034564433310
0.010279663295576
Iteration 3
0.000000000
1.037242198
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
x
0 = 5.227918429
λ
0 = 0.036588313752304
Iteration 1
0.143632267
1.025322691
5.087833489750430
0.038630846139210
Iteration 2
−0.000007263
1.025426533
5.087840573092400
0.038630738574793
Iteration 3
0.000000000
1.025426528
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
5.087840573092420
0.038630738574792Table 5
Halley procedure. Option 5: fixed initial starting pointx0 = 7.990256504→λ0 = 0.015663210285978339 from Section 2.2.2, indirect calculation of λ through the transmission factor x: Equation (13), and the symbolic derivatives f′(x) and f″(x): Equations (11) and (14).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
−1.945484250
1.040493788
−0.001887828
9.863203600915390
0.010279310950983
Iteration 2
0.000175332
1.037241928
−0.001596798
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
2.973246188
1.023435376
−0.000632309
5.087698791122220
0.038632891696967
Iteration 2
−0.000145387
1.025426633
−0.000744326
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
5.087840573092420
0.038630738574792Table 6
Third-order Householder’s procedure. Option 6: starting pointx0 depends on input parameters: Equation (2), indirect calculation of λ through the transmission factor x: Equation (15), and the symbolic derivatives f′(x), f″(x), and f‴(x): Equations (11), (14), and (16).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
f‴(x), Equation (16)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.495092014
1.036495031
−0.001533392
0.000128855
9.863034531578420
0.010279663364062
Iteration 2
−0.000000034
1.037242198
−0.001596821
0.000136933
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
0.000136933
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
f‴(x), Equation (16)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.143632267
1.025322691
−0.000738253
0.000043046
5.087840573035260
0.038630738575660
Iteration 2
0.000000000
1.025426528
−0.000744320
0.000043578
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
0.000043578
5.087840573092420
0.038630738574792Table 7
Schröder procedure. Option 7: fixed initial starting pointx0 = 7.990256504→λ0 = 0.015663210285978339 from Section 2.2.2, indirect calculation of λ through the transmission factor x: Equation (17), and the symbolic derivatives f′(x) and f″(x): Equations (11) and (14).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
−1.945484250
1.040493788
−0.001887828
9.863198212166060
0.010279322183170
Iteration 2
0.000169742
1.037241937
−0.001596799
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
2.973246188
1.023435376
−0.000632309
5.087701128882780
0.038632856193927
Iteration 2
−0.000142990
1.025426632
−0.000744326
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
5.087840573092420
0.038630738574792


## 3.2.1. Indirect Calculation ofλ through the Transmission Factor x=1/λ with the Derivative Calculated Analytically

The first derivativef′ with respect to x can be obtained analytically, Equation (9) (also Equation (11) gives the same results):(9)f′x=ddxfx=1+2⋅2.51/Re⋅ln10/ε/3.7⋅D+2.51/Re⋅x︸1st  derivativef′x−analytical.Initial value of the flow friction factorλ0 should be chosen and the residue f/f′ is calculated in order to start the Newton–Raphson procedure:(10)x1=x0−fx0/f′x0︸Newton−Raphson=x0−x0+2⋅log102.51⋅x0/Re+ε/3.7⋅D1+2⋅2.51/Re⋅ln10/ε/3.7⋅D+2.51/Re⋅x0.The procedurexi+1 = xi − f(xi)/f′(xi) should be followed until the residue f(xi)/f′(xi) ≈ 0. Then, the final solution is λn=xn−2, 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 x0 = 6.445695939 for Re = 5·106, ε/D = 2.5·10−5).Table 3
Newton–Raphson procedure. Option 3: fixed initial starting pointx0 = 6.445695939→λ0 = 0.024069128765101 from Section 2.2.1, indirect calculation of λ through the transmission factor x: Equation (10), and analytical derivative f′(x): Equation (9).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (9)
x
0 = 6.445695939
λ
0 = 0.024069128768719
Iteration 1
−3.554956085
1.043635910
9.852014225862620
0.010302673560706
Iteration 2
−0.011430857
1.037259804
9.863034470914730
0.010279663490514
Iteration 3
−0.000000097
1.037242198
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (9)
x
0 = 6.445695939
λ
0 = 0.024069128768719
Iteration 1
1.391712393
1.024454486
5.087204750239650
0.038640395682209
Iteration 2
−0.000651990
1.025427001
5.087840572945700
0.038630738577020
Iteration 3
0.000000000
1.025426528
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.046830475
5.087840573092420
0.038630738574792


## 3.2.2. Indirect Calculation ofλ through the Transmission Factor x=1/λ with the Symbolic Derivative

The exact analytical expression of the first derivativef′ with respect to x can be obtained in MATLAB, Equation (11); results are the same as using Equation (9):(11)f′x=ddxfx=5.02/Re⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re+1=9287⋅ln10⋅x+1000⋅ln10⋅ε/D⋅Re+18574/ln10⋅9287⋅x+1000⋅ε/D⋅Re︸1st  derivativef′x−MATLAB.Initial value of the flow friction factorλ0 should be chosen and the residue f/f′ is calculated in order to start the Newton–Raphson procedure:(12)x1=x0−fx0/f′x0︸Newton−Raphson=x0−x0+2⋅log102.51⋅x0/Re+ε/3.7⋅D21384.11⋅x0+2302.58⋅ε/D⋅Re+18574/21384.11⋅x0+2302.58⋅ε/D⋅Re.The procedurexi+1 = xi − f(xi)/f′(xi) should be followed until the residue f(xi)/f′(xi) ≈ 0. Then, the final solution is λn=xn−2, where n = i + 1 is the final iteration.The iterative procedure can be accelerated using Halley’s formula instead of the Newton–Raphson:(13)x1=x0−fx0/f′x0/1−fx0·f″x0/2⋅f′x02︸Halley=x0−2⋅fx0⋅f′x0/2⋅f′x02−fx0⋅f″x0︸Halley.In general,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is the final iteration in which xn ≈ xn+1.The second derivativef″(x) in respect of x is required:(14)f″x=ddxf′x=−12.6/Re2⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re2=−172496738/ln10⋅9287⋅x+1000⋅ε/D⋅Re2︸2ndderivativef″x−MATLAB.The Newton–Raphson method belongs to the 1st order and the Halley to the 2nd order of Householder’s method, while the 3rd order can be expressed using the following equation:(15)x1=x0−6⋅fx0⋅f′x02−3⋅fx02⋅f″x0/6⋅f′x03−6⋅fx0⋅f′x0⋅f″x0+fx02⋅f‴x0︸3rdorderHouseholder.Again,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is the final iteration in which xn ≈ xn+1.The required 3rd derivative f‴(x) can be expressed using the following equation:(16)f‴x=ddxf″x=63.253/Re3⋅ln10⋅10/37⋅ε/D+251⋅x/100⋅Re3=32039544116.2/ln10⋅9287⋅x+1000⋅ε/D⋅Re3︸3rd  derivativef‴x−MATLAB.Also here one has to be underlined that the Halley method [74] is not the unique Householder method of the 2nd order [22]. For example, the Schröder method [25] also belongs to the group:(17)x1=x0−fx0/f′x0−f″x0⋅fx02/2⋅f′x03︸Schröder.Further,x1 = xi+1 and x0 = xi; i = 0 to n, where n + 1 is final iteration in which xn ≈ xn+1.Using the presented Householder procedures, the 1st order: the Newton–Raphson, the 2nd order: Halley, and the 3rd order, the unknown flow friction factor λ should be calculated for the two given pairs of the Reynolds number Re and the relative roughness ε/D: (1) (Re = 5·106, ε/D = 2.5·10−5)⟶λ = 0.010279663295529 and (2) (Re = 3·104, ε/D = 9·10−3)⟶λ = 0.038630738574792. The calculation presented in Tables 4–7 is through the transmission factor x, using the symbolic derivative f′(x), but with the different initial starting point λ0.Table 4
Newton–Raphson procedure. Option 4: starting pointx0 depends on input parameters: Equation (2), indirect calculation of λ through the transmission factor x: Equation (12), and the symbolic derivative f′(x): Equation (11).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.495092014
1.036495031
9.862863625818000
0.010280019623455
Iteration 2
−0.000177305
1.037242471
9.863034564433310
0.010279663295576
Iteration 3
0.000000000
1.037242198
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
x
0 = 5.227918429
λ
0 = 0.036588313752304
Iteration 1
0.143632267
1.025322691
5.087833489750430
0.038630846139210
Iteration 2
−0.000007263
1.025426533
5.087840573092400
0.038630738574793
Iteration 3
0.000000000
1.025426528
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
5.087840573092420
0.038630738574792Table 5
Halley procedure. Option 5: fixed initial starting pointx0 = 7.990256504→λ0 = 0.015663210285978339 from Section 2.2.2, indirect calculation of λ through the transmission factor x: Equation (13), and the symbolic derivatives f′(x) and f″(x): Equations (11) and (14).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
−1.945484250
1.040493788
−0.001887828
9.863203600915390
0.010279310950983
Iteration 2
0.000175332
1.037241928
−0.001596798
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
2.973246188
1.023435376
−0.000632309
5.087698791122220
0.038632891696967
Iteration 2
−0.000145387
1.025426633
−0.000744326
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
5.087840573092420
0.038630738574792Table 6
Third-order Householder’s procedure. Option 6: starting pointx0 depends on input parameters: Equation (2), indirect calculation of λ through the transmission factor x: Equation (15), and the symbolic derivatives f′(x), f″(x), and f‴(x): Equations (11), (14), and (16).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
f‴(x), Equation (16)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.495092014
1.036495031
−0.001533392
0.000128855
9.863034531578420
0.010279663364062
Iteration 2
−0.000000034
1.037242198
−0.001596821
0.000136933
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
0.000136933
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
f‴(x), Equation (16)
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
0.143632267
1.025322691
−0.000738253
0.000043046
5.087840573035260
0.038630738575660
Iteration 2
0.000000000
1.025426528
−0.000744320
0.000043578
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
0.000043578
5.087840573092420
0.038630738574792Table 7
Schröder procedure. Option 7: fixed initial starting pointx0 = 7.990256504→λ0 = 0.015663210285978339 from Section 2.2.2, indirect calculation of λ through the transmission factor x: Equation (17), and the symbolic derivatives f′(x) and f″(x): Equations (11) and (14).
Re = 5·106, ε/D = 2.5·10−5
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
−1.945484250
1.040493788
−0.001887828
9.863198212166060
0.010279322183170
Iteration 2
0.000169742
1.037241937
−0.001596799
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
1.037242198
−0.001596821
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(x), Equation (8)
f′(x), Equation (11)
f″(x), Equation (14)
x
0 = 7.990256504
λ
0 = 0.015663210285978
Iteration 1
2.973246188
1.023435376
−0.000632309
5.087701128882780
0.038632856193927
Iteration 2
−0.000142990
1.025426632
−0.000744326
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000000
1.025426528
−0.000744320
5.087840573092420
0.038630738574792


## 3.3. Secant Method

The secant method is similar to the Newton–Raphson; it requires two starting pointsλ0 and λ−1 but does not require calculation of the derivatives [71]. The approach with the direct calculation of λ with the two required starting points λ0 and λ−1 is formulated:(18)λ1=λ0−fλ0fλ−1−fλ0/λ−1−λ0.Counteri starts from i − 1 and goes to n + 1 in which λn = λn+1.The approach through the transmission factor x is as follows:(19)x1=x0−fx0fx−1−fx0/x−1−x0.As already described, counteri also starts from i − 1 and goes to n + 1 in which xn = xn+1.The flow friction factorλ is calculated in Tables 8 and 9 for two pairs of the Reynolds number and the relative roughness (1) Re = 5·106, ε/D = 2.5·10−5, and (2) Re = 3·104, ε/D = 9·10−3, using the secant procedure with direct calculation of λ and indirect through the transmission factor x.Table 8
Secant procedure. Option 8: two initial starting pointsλ0 and λ−1 required: starting point λ−1 is with fixed value λ−1 = 0.024069128765101 (i.e., x−1 = 6.445695939) as in Section 2.2.1, while starting point λ0 depends on input parameters: Equation (2), and direct calculation of λ: Equation (18).
Re = 5·106, ε/D = 2.5·10−5
f(λi−1), Equation (5)
f(λi), Equation (5)
f
λ
i
−
1
−
f
λ
i
/
λ
i
−
1
−
λ
i
λ
−1 = 0.024069128765101
λ
0 = 0.009352225155363
Iteration 1
0.495092014
−3.554956084
−275.1970255
0.011151270814558
Iteration 2
−0.408071981
0.495092014
−502.0239429
0.010338417191085
Iteration 3
−0.029111936
−0.408071981
−466.2094551
0.010275973292109
Iteration 4
0.001836644
−0.029111936
−495.6221591
0.010279679026163
Iteration 5
−0.000007828
0.001836644
−497.734448
0.010279663299743
Iteration 6
−0.000000002
−0.000007828
−497.6011214
λ = 0.010279663295529
Control step
0.000000000
−0.000000002
−497.6012179
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(λi−1), Equation (5)
f(λi), Equation (5)
f
λ
i
−
1
−
f
λ
i
/
λ
i
−
1
−
λ
i
λ
−1 = 0.024069128765101
λ
0 = 0.036588313752304
Iteration 1
1.391712394
0.143632267
−99.69340079
0.038029053721052
Iteration 2
0.143632267
0.041110009
−71.15944585
0.038606770549177
Iteration 3
0.041110009
0.001619232
−68.35663251
0.038630458556837
Iteration 4
0.001619232
0.000018909
−67.5583902
0.038630738444645
Iteration 5
0.000018909
0.000000009
−67.52699052
λ = 0.038630738574792
Control step
0.000000009
0.000000000
−67.52662212
0.038630738574792Table 9
Secant procedure. Option 9: two initial starting pointsx0 and x−1 required: starting point λ−1 is with fixed value x−1 = 6.445695939 (i.e., λ−1 = 0.024069128765101) as in Section 2.2.1, while starting point λ0 depends on input parameters: Equation (2), and indirect calculation of λ through the transmission factor x: Equation (19).
Re = 5·106, ε/D = 2.5·10−5
f(xi−1), Equation (8)
f(xi), Equation (8)
f
x
i
−
1
−
f
x
i
/
x
i
−
1
−
x
i
x
−1 = 6.445695939
λ
−1 = 0.024069128765101
x
0 = 10.34052343
λ
0 = 0.009352225155363
Iteration 1
−3.554956084
0.495092014
1.039853012
9.864406125318800
0.010276804896656
Iteration 2
0.495092014
0.001422639
1.03686501
9.863034066961850
0.010279664332547
Iteration 3
0.001422639
−0.000000516
1.037241104
9.863034564456330
0.010279663295528
Iteration 4
−0.000000516
0.000000000
1.0372422
x = 9.863034564455800
λ = 0.010279663295529
Control step
0.000000000
0.000000000
1.037162162
9.863034564455800
0.010279663295529
Re = 3·104, ε/D = 9·10−3
f(xi−1), Equation (8)
f(xi), Equation (8)
f
x
i
−
1
−
f
x
i
/
x
i
−
1
−
x
i
x
−1 = 6.445695939
λ
−1 = 0.024069128765101
x
0 = 5.227918429
λ
0 = 0.036588313752304
Iteration 1
1.391712394
0.143632267
1.024883541
5.087773465040530
0.038631757665255
Iteration 2
0.143632267
−0.000068814
1.025374564
5.087840576494990
0.038630738523123
Iteration 3
−0.000068814
0.000000003
1.025426553
x = 5.087840573092420
λ = 0.038630738574792
Control step
0.000000003
0.000000000
1.025426591
5.087840573092420
0.038630738574792Calculation 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·106, ε/D = 2.5·10−5, and the required number of iterations is 4 in direct calculation and 3 in indirect for Re = 3·104, ε/D = 9·10−3.


## 3.4. Three-Point Methods

Three-point iterative methods require in every iteration an evaluation of the functionf at three points: x0, y0, and z0. However, these methods converge very fast to the accurate solution. The mathematical background of the three-point iterative methods is given by Sharma and Arora [28]. Here, we will apply the Džunić–Petković–Petković three-point iterative method to the Colebrook equation [26, 28]. It requires only one iteration (up to 2 in the worst cases) to reach the final accurate solution and we will show all steps to calculate the friction factor λ, Equation (20), where the numerical values are for Re = 5·106, ε/D = 2.5·10−5. The initial starting point is x0 = 7.273124147 as described in Section 2.2.3.(20)x0=7.273124147,fx0=x0+2⋅log102.51⋅x0Re+ε3.7⋅D=−2.692152546,f′x0=9287⋅ln10⋅x0+1000⋅ln10⋅ε/D⋅Re+18574ln10⋅9287⋅x0+1000⋅ε/D⋅Re=5.02Re⋅ln10⋅10/37⋅ε/D+2.51⋅x0/Re+1=1.041894438,y0=x0−fx0f′x0=7.273124147−−2.6921525461.041894=9.857025593360860,fy0=y0+2⋅log102.51⋅y0Re+ε3.7⋅D=−0.006232787,z0=y0−fx0fx0−2·fy0·fy0f′x0=9.863035589,fz0=z0+2⋅log102.51⋅z0Re+ε3.7⋅D=−0.006232787,x1=z0−fz0f′x0·1−2·fy0/fx0−fy0/fx02·1−fz0/fy0·1−2·fz0/fx0=9.863034564,x1=9.863034564⟶λ1=0.010279663295529.Three-point methods are also shown in detail in Praks and Brkić [75].


## 3.5. Expressed through the Lambert W-Function

The mathematical background of the Lambert W-function is presented in the literature [44–46]. Moreover, an application of the Lambert W-function in hydraulic problems is given by [41, 54, 73].The Lambert W-functionW(y) [46] is the solution of y=z⋅ez, which needs to be in the appropriate form:(21)fz=z⋅ez−y=0.The first derivativef′(z) is as follows:(22)f′z=ez⋅z+1.Choose initial valuez0, calculate the residue f/f′, and start the procedure:(23)z1=z0−fz0f′z0=z0−z0⋅ez0−y/ez0⋅z0+1︸Newton−Raphson=z0−z0⋅ez0−Re/2.18ez0⋅z0+1.Then follow the procedurezi+1 = zi − f(zi)/f′(zi) until the residue f(zi)/f′(zi)≈0, where n = i + 1 in final iteration.Halley’s procedure:(24)z1=z0−fz0f′z0−fz0⋅f″z0/2⋅f′z0=z0−z0⋅ez0−y/ez0⋅z0+1−z0⋅ez0−y⋅z0+2/2⋅z0+1︸Halley,where the second derivative is(25)f′z=ez⋅z+2.The Schröder expression is as follows:(26)z1=z0−fz0f′z0−f″z0⋅fz022⋅f′z03=z0−z⋅ez−y/ez⋅z+1−ez⋅z+2⋅z⋅ez−y2/2⋅ez⋅z+13︸Schröder.Further, in all cases, the Newton–Raphson, the Halley, and the Schröder;z1 = zi+1 and z0 = zi; i = 0 to n, where n + 1 is final iteration in which zn ≈ zn+1.The argument of the Lambert W-functiony, in our case defined by Equation (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 the Halley method for Re = 5·106 and Re = 3·104 where initial starting point is set as z0 = 15 as recommended in Section 2.3 of this paper (for z0 < 8.814, the Newton–Raphson procedure cannot start).Table 10
Calculation of W(y) wherey = Re/2.18 using the Newton–Raphson, the Halley, and the Schröder iterative methods.
Re = 5·106, y = Re/2.18 = 2293411.45
Re = 3·104, y = Re/2.18 = 13760.47
Newton–Raphson
Halley
Schröder
Newton–Raphson
Halley
Schröder
Iteration 0
z
0 = 15
z
0 = 15
z
0 = 15
z
0 = 15
z
0 = 15
z
0 = 15
Iteration 1
14.10634749
13.29860556
13.68208338
14.06276308
13.1333396
13.59610616
Iteration 2
13.28604947
12.2757343
12.62556802
13.12986539
11.29171838
12.20340124
Iteration 3
12.62863905
12.14855784
12.17738057
12.20257063
9.520829163
10.83006437
Iteration 4
12.25343232
z = 12.14835704
12.14836628
11.28354302
8.068323472
9.505729616
Iteration 5
12.15407754
12.14835704
z = 12.14835704
10.37904335
7.530266826
8.341483562
Iteration 6
12.14837461
—
12.14835704
9.504505014
7.512930233
7.637280252
Iteration 7
z = 12.14835704
—
—
8.697314341
z = 7.512929679
7.513654122
Iteration 8
12.14835704
—
—
8.037456295
7.512929679
z = 7.512929679
Iteration 9
—
—
—
7.640105762
—
7.512929679
Iteration 10
—
—
—
7.52154464
—
—
Iteration 11
—
—
—
7.512971011
—
—
Iteration 12
—
—
—
7.51292968
—
—
Iteration 13
—
—
—
z = 7.512929679
—
—
Iteration 14
—
—
—
7.512929679
—
—


## 4. Approximations: Simplified Equations for Engineering Practice

Using the optimal fixed initial starting point for the Halley and the Schröder method as explained in Section2.2.2, the first iteration of the procedures from Section 3.2.2, the simplification using the fact that the first derivative of the Colebrook function is almost always near one, f′ ≈ 1, and using acceleration through Equation (28) [76, 77], the following approximations, Equation (27), can be formed. Using Equation (27), the maximal relative error in the domain of applicability of the Colebrook equation is 8.29% (Figure 6), and using acceleration Equation (28), that is, single fixed-point iterative method [12], the maximal relative error is 0.69% (Figure 7), 0.0617% (Figure 8), etc.(27)1λ0≈8−2⋅A/2−A·B︸Halley≈8−A−A2⋅B/2︸Schröder≈8−6⋅A−3⋅A2⋅B/6−6⋅A⋅B+A2⋅C︸3rd  order.Figure 6
Distribution of the relative error over the domain of applicability of the Colebrook equation of the approximation, Equation (27); for i = 0, the maximal relative error is 8.29%.Figure 7
Distribution of the relative error over the domain of applicability of the Colebrook equation of the approximation, Equation (28); for i = 0, the maximal relative error is 0.69%.Figure 8
Distribution of the relative error over the domain of applicability of the Colebrook equation of the approximation, Equation (28); for i = 1, the maximal relative error is 0.0617%.Then,λ0 from Equation (27) is used in(28)1λ1≈−2⋅log102.51/Re⋅1/λ0+ε/3.7⋅D︸1stColebrook'sacceleration,1λ2≈−2⋅log102.51/Re⋅1/λ1+ε/3.7⋅D︸2stColebrook'sacceleration=−2⋅log102.51/Re⋅−2⋅log102.51/Re⋅1/λ0+ε/3.7⋅D︸1stColebrook'sacceleration+ε/3.7⋅D︸2stColebrook'sacceleration,⋮1λi+1≈−2⋅log102.51/Re⋅1/λi+ε/3.7⋅D︸Colebrook'sacceleration,whereA, B, and C are(29)A≈8+2⋅log1016Re+ε3.7⋅D,B≈−74914381.46∇2,C≈1391459721232.67∇3,where ∇(30)∇≈74205.5+1000⋅εD⋅Re.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% can be reached using only three logarithmic forms (in all cases, exponents are only whole numbers) [6, 8, 33].The error analysis is in the domain of applicability of the Colebrook equation and can be further reduced using one more accelerating step as shown, genetic algorithms [31, 78, 79], Excel fitting tool [80], or Monte Carlo [81, 82].Zigrang and Sylvester [77] used the similar approach. They use an iterative procedure to produce an approximation, Equation (31), where α1 and α2 represent internal iterative steps, while constant 13 in α2 is fixed starting point. Shacham [83] used the same approach, but he set the fixed starting point as 14.5, while [31] also optimized that parameter which makes clear that it does not to be fixed.(31)1λ≈−2·log10ε3.7·D−5.02Re·α1,α1≈log10ε3.7·D−5.02Re·α2,α2≈log10ε3.7·D−13Re.Equation (31) includes three logarithmic forms and produces the maximal relative error of 0.114%. The same approach is used by Schorle et al. [83]. However, although Equation (31) and also the here presented Equation (28) include three logarithmic forms, the maximal relative error of Equation (28) for i = 1 is only 0.0617%. Thus, the here presented approach reduced the maximal relative error of the three logarithmic form approximation by factor 0.114/0.0617∼1.8.


## 5. Conclusions

The paper presents a fast but reliable approximations and iterative methods for pipeline hydraulics, useful for the reliable modelling of water and gas distribution networks where a large number of network simulations of random component failures and their combinations need to be automatically evaluated and statistically analyzed [84–89]. Accurate, fast, and reliable estimation of the flow friction factor is essential for the evaluation of pressure drops and flows in large network of pipes, because, for example, compression station failures can be only approximated in the transmission level by user-defined logic rules obtained from hydraulic software [81, 82, 90–92]. Iterative solutions and approximations for the calculation of the flow friction factor are implemented in software packages which are in common use in everyday engineering practice [88]. So in this paper, we analyzed selected iterative procedures in order to solve the Colebrook equation [93, 94], and we found that up 2 to 3 iterations of the Halley and the Schröder method are suitable for the accuracy required by engineering practice, when the fixed initial starting point described in Section 2.2.2 is applied. On the other hand, using a three-point iterative method with the same initial conditions, the required high accuracy can be reached after only 1 iteration (2 in the worst case) but using three internal steps [26–28].Moreover, for implementing simplified engineering calculations, we recommend taking into consideration the following:(1)
The indirect calculation ofλ through the transmission factor x in most cases accelerates the iterative procedures.(2)
Knowing that the Colebrook equation is used in engineering practice only in the limited domain of the Reynolds number Re between 4000 and 108, and for relative roughness of inner pipe surface ε/D up to 0.05, we evaluated the number of iterations required to reach sufficient accuracy. We detected zones of input parameters, in which iterative methods converge slowly, and thus, additional number of iterations is required. Therefore, we put the fixed initial point to start the iterative procedure in those zones in order to decrease the number of required iterations.(3)
Using the simplified Halley and the simplified Schröder procedure with the fixed starting point, after only one iteration, one can reach results with good accuracy (error up to 8.29%). This is near the accurate value, and therefore, the simplified Newton–Raphson method can be used from the second iteration to reach an accuracy of 0.69%. Moreover, one can reach an accuracy of 0.0617% using the third iteration. With the methods proposed herein, the first derivative of the Colebrook function is always near one,f′≈1 (for f′⟶1, the Newton–Raphson method becomes the fixed-point method). Consequently, accuracy of 0.69% can be reached using only two logarithmic forms and of 0.0617% with only three logarithmic forms without requiring extensive computational efforts (the goal being to use the least possible number of logarithmic functions or functions with noninteger power) [6, 33]. Moreover, the computational cost of iterations can also be reduced by Padé polynomials [95, 96]. We also analyzed the Colebrook equation expressed through the Lambert W-function, and we found that the Halley and the Schröder method can be advised in comparison with the Newton–Raphson method (where the problem with the initial starting point exists).The Colebrook equation is valid only for turbulent flow, while to include also laminar regime, different models should be used [98].



---

*Source: 5451034-2018-12-10.xml*