What Can Students Learn While Solving Colebrook’s Flow Friction Equation?

: Even a relatively simple equation such as Colebrook’s o ﬀ ers a lot of possibilities to students to increase their computational skills. The Colebrook’s equation is implicit in the ﬂow friction factor and, therefore, it needs to be solved iteratively or using explicit approximations, which need to be developed using di ﬀ erent approaches. Various procedures can be used for iterative methods, such as single the ﬁxed-point iterative method, Newton–Raphson, and other types of multi-point iterative methods, iterative methods in a combination with Pad é polynomials, special functions such as Lambert W, artiﬁcial intelligence such as neural networks, etc. In addition, to develop explicit approximations or to improve their accuracy, regression analysis, genetic algorithms, and curve ﬁtting techniques can be used too. In this learning numerical exercise, a few numerical examples will be shown along with the explanation of the estimated pedagogical impact for university students. Students can see what the di ﬀ erence is between the classical vs. ﬂoating-point algebra used in computers.


Introduction
The Colebrook equation for flow friction is an empirical formula developed by C. F. Colebrook in 1939 [1] and it is based on his joint experiment, which was performed with Prof. C. M. White in 1937 [2]. Colebrook was at this time a PhD student of Prof. White at Imperial College in London, UK [3]. The experiment was performed with a set of pipes of which the inner side was covered with sand of different grain sizes, while one was left completely without sand to be smooth. The sand was fixed on the pipe's inner surface using a sort of glue as adhesive material, while air was used as the fluid in this experiment.
The Colebrook equation follows logarithmic law and, at first sight, it seems to be very simple. However, it is given implicitly, with respect to the unknown flow friction factor f in a way that it can be expressed explicitly only in terms of the Lambert W function, with further difficulties in its evaluation [4][5][6] (it can be rearranged for gas flow [7]). Otherwise, it needs to be solved iteratively [8][9][10] or explicitly, using approximate formulas specially developed for such a purpose [10][11][12][13][14][15][16][17]. The unknown flow friction factor f in the Colebrook's equation depends on the known Reynolds number Re and the relative roughness of the inner pipe surface ε D , all three of which are dimensionless, shown in Equation (1): 1

2.51
Re (2) Re = 4.6 × 10 7 and ε D = 0.037. Corresponding solutions for the flow friction factor are (1) f = 0.016050961, and (2) f = 0.062427396, respectively. Positions of these two examples with respect to the hydraulic regimes (laminar, transition to turbulent, and turbulent, which includes smooth, smooth to rough, and fully turbulent rough regime) are shown in Figure 1. Cases that can capture the nonlinearity of the Colebrook equation are: (a) Smooth pipe at very low, moderate, and very high Re, (b) pipe with moderate roughness at very low, moderate, and very high Re, and (c) pipe with high roughness ( ε D ∼ 0.05) and again, low, moderate, and high Re. Example 1 is in the zone of moderate values of the Reynolds number Re, where the proposed value of relative roughness, which is small ( ε D ∼ 10 −4 ), is large enough to provide a transition from smooth towards the turbulent regime. On the other hand, Example 2 with higher values of Re and ε D is in a fully turbulent zone. Using Figure 1, other different examples can be provided for these teaching exercises.
Example 1 is set to be exactly in the zone where the smooth nature of the viscose zone of the inner pipe surface starts to be too thin and when roughness starts to have influence [47,48]. Example 1 is set to be exactly in the zone where the smooth nature of the viscose zone of the inner pipe surface starts to be too thin and when roughness starts to have influence [47,48]. From Figure 1, students should see that our numerical Example 1 is set to be in the turbulent zone. More precisely, Example 1 is set at the end of the smooth turbulent regime, where the transition to fully turbulent begins. On the contrary, Example 2 is located in the fully developed rough turbulent zone. Students should understand the physical interpretation of different hydraulic regimes.
For this task, 45 min should be assigned.

Trial/Error Method
To do this task, students should guess somehow the initial value of the friction factor , then insert it into the Colebrook equation, and finally, calculate a new value of until an error criterion is fulfilled. Using this iterative procedure, students should show how many iterations will probably be required, as the number of iterations depends on the initial guess and whether such a rudimentary iterative procedure will always lead to the final balanced solution or not. Here will be introduced a concept of convergence, divergence, and oscillations in iterative procedures. Students should understand how significant it is to find the adequate value for the starting point for an iterative procedure [9]. For our Example 1, which uses = 2.3 × 10 and = 10 , the starting point needs to be chosen between −2.469 and 1574003 when Microsoft (MS) Excel is used and students should discuss why [10,20]. Subsequently, students should analyze the starting point for Example 2, in which = 4.6 × 10 and = 0.037. Students should explain consequences of the inappropriately chosen initial starting point. Finally, students should understand that the starting From Figure 1, students should see that our numerical Example 1 is set to be in the turbulent zone. More precisely, Example 1 is set at the end of the smooth turbulent regime, where the transition to fully turbulent begins. On the contrary, Example 2 is located in the fully developed rough turbulent zone. Students should understand the physical interpretation of different hydraulic regimes.
For this task, 45 min should be assigned.

Trial/Error Method
To do this task, students should guess somehow the initial value of the friction factor f , then insert it into the Colebrook equation, and finally, calculate a new value of f until an error criterion is fulfilled. Using this iterative procedure, students should show how many iterations will probably be required, as the number of iterations depends on the initial guess and whether such a rudimentary iterative procedure will always lead to the final balanced solution or not. Here will be introduced a concept of convergence, divergence, and oscillations in iterative procedures. Students should understand how significant it is to find the adequate value for the starting point for an iterative procedure [9]. For our Example 1, which uses Re = 2.3 × 10 5 and ε D = 10 −4 , the starting point needs to be chosen between −2.469 and 1574003 when Microsoft (MS) Excel is used and students should discuss why [10,20]. Subsequently, students should analyze the starting point for Example 2, in which Re = 4.6 × 10 7 and ε D = 0.037. Students should explain consequences of the inappropriately chosen initial starting point. Finally, students should understand that the starting point for −2.5, for example, will actually be 3.71 − 2.5 = −2.499973046 for our Example 1, in which ε D = 10 −4 (the negative starting point does not have physical meaning and it is chosen purely for testing methods and to introduce the concept of "thinking outside of the box", which could contribute to development of creative thinking in general.
For this task, 45 min should be assigned.

Fixed-Point Iterative Procedure
The fixed-point iterative method is suitable for spreadsheet solvers, such as MS Excel [10,20]. To prepare this task, students should rearrange the Colebrook equation in the suitable form F(x) = 0; The iterative process starts with a selected starting point x 0 and continues by x i+1 = x i − F(x i ) until a stopping criteria is fulfilled, for example, when the results from two subsequent iterations are almost equal ( x i+1 − x i+1 ∼ 0). From this approach, students will learn that the trial/error method can be used in a more systematic way. This task can improve the ability of students to use spreadsheet solvers, because students should set the solver to allow iterative computations. Moreover, the maximal number of iterations and the stopping criteria represented by the minimal value between two successive iterations should be set in the solver. For example, the MS Excel spreadsheet solver allows a maximal number of iterations of 32,767, and students should figure out why this constraint is set. The answer is that 2 15 = 32,768 as this represents the maximum it can be set for 16-bit registers in binary logic because one bit is reserved for the sign (±). In that way, students will learn that for human perception, more common is the decimal system, which is used, for example, for banknotes, $1, $10, $100, while, for computers, the binary system is more suitable, which commonly uses 2, 4, 8, 16, 32, 64, etc.
In our case, (1) Re = 2.3 × 10 5 and ε .71 + 1,574,003 ≈ 1,574,003, the number of required iterations for nine decimal accuracy will be 11 in both cases, i.e., the accurate solution will be 1 √ f 11 = 7.8931. Students should make further tests with the same example, but with different starting points and, in addition, perform a similar test for Example 2.
Using the spreadsheet solver, students should first solve the Colebrook equation iteratively, using only one cell of the Colebrook equation for smooth conditions, when ε D is zero. Students should realize how important the role of the initial starting point is, because, in this case, software will declare a division by zero error. In other conditions, when ε D 0, the initial starting point will be automatically set as f 0 = log 10 ε D

.
Here, students may try to make a sketch of the function, i.e., to make a two-dimensional (2D) diagram of the Colebrook equation. Here, two possibilities should be analyzed: 1 as suggested by Rose [46], or f = F Re, ε D as suggested by Moody [18] (approach from Figure 1 of this article can be used). In addition, students can analyze whether a simple transformation x = 1 √ f accelerates solving of the Colebrook equation. For this task, 45 min should be assigned.

Newton-Raphson Iterative Procedure
In this part, students will learn how to find the first derivative of the function with respect to the variable flow friction factor f , where x = 1 √ f and F(x) = x + 2 log 10

2.51
Re ·x + ε D 3.71 , as in Equation (2): Re· ln(10)· 2.51 The procedure for finding the first derivative can be done manually, but also for Equation (2), it can be generated using appropriate software, e.g., Matlab by the "diff" command.
The first derivative is used in the Newton-Raphson procedure as and students may test whether the Newton-Raphson procedure is faster or not, compared with the simple fixed-point 3.71 . Students should realize that the previous example with the fixed point method is set for solving 3.71 + 1,574,003 ≈ 1,574,003 → f 0 = 4.04 × 10 −13 , which will give accurate results. Unfortunatelly, this iteration process 173611 will report a division by zero error in the fourteenth iteration and the calculation will stop without obtaining the precise solution. Here, students should understand that even a more powerful iterative procedure can fail when the problem is not set in an appropriate way and when the role of a good initial starting point is not fully understood.
The Newton-Raphson method, in which F ( f i ) = 1 is assumed, is equivalent to the fixed-point method in the form f i+1 = f i − F( f i ). However, this formulation shows very bad convergence properties. Students should understand why the iterative process while the iterative process f i+1 = f i − F( f i ) is slow. Students should again analyze coordinates in diagrams from Rose and Moody, and discuss which formulation is better [18].
Here, it should be noted that the function F(x) = x + 2 log 10

3.71
in form shows very slow convergence when it is used in the iterative Students should discover why. Students should investigate how the initial starting point affects the required number of iterations.
For this task, 90 min should be assigned (including repetition and clarifications from the previous tasks).

Multi-Point Iterative Procedures
The shown Newton-Raphson method belongs to the two-point iterative methods and, in this task, students should find more powerful iterative methods. For example, the proposed methods for this exercise are three-point iterative methods, such as Halley or Schröder or even more powerful multi-point methods [8,9]. Students should figure out whether the number of required iterations to find the accurate solution decreases with the increased complexity of the chosen method. If the number of required iterations does not decrease, students should analyze why. In addition, it needs to be found how important it is to find the adequate initial starting point. In this task, students should find the second derivative of the Colebrook equation, which will increase their computing skills.
In this exercise, students should develop the ability to make a balance between complexity, speed, and required accuracy. Students obtain an experience for choosing the appropriate iterative method.
Here, we will not repeat the Newton-Raphson calculations in the form of Equation (3), which gives fast accurate results, but we will show the Halley method: Equation (4) for which the first derivative F (x) and the second F (x), with respect to variable x = 1 √ f i , are given by Equations (5) and (6), respectively.
For the first derivatives, [8,9,[21][22][23] should be consulted. For Example 1, in which Re = 2.3 × 10 5 and ε D = 10 −4 , for initial starting point .71 + 1574003 ≈ 1574003, the final balanced solution will be achieved after 27 iterations; .89313398624169 → f 27 = 0.016050961 . However, for the initial starting point , the procedure in MS Excel will report a division by zero error in the third step, so it will be a good opportunity for students to reexamine the importance of the appropriate chosen starting point. Thus, x 0 = 1 √ f 0 , which in this case can be selected from −1.172 to 1,583,778.
Of course, negative starting points do not represent physical meaning. A similar exercise should be done for Example 2, in which Re = 4.6 × 10 7 and ε D = 0.037. For this task, at least 90 min should be assigned, as multi-point iterative methods are more complex than two-point methods.

Special Iterative Methods
While simple operations, such as adding, subtracting, multiplication, and division, require minimal system resources of computers, more complex operations, such as logarithmic, non-integer power, and similar function calls, require multiple floating point operations to be executed in the Central Processing Unit (CPU), which represents a time-and energy-consuming task [15,[25][26][27][28]. For example, students should understand that the execution of one non-integer power operation in a computer goes through execution of one logarithmic and one exponential function. Some numbers are too big or too small to be placed in registers of computers, which have 32 or 64 bits [37,38].

Padé Approximant
As noted, one iterative cycle involves the evaluation of one logarithmic function when solving the Colebrook equation. The Colebrook's law is in its nature logarithmic, but here, students will learn how to avoid multiple evaluations of the logarithmic function by replacing it with a Padé approximant in all iterations, with the exception of the first one [30]. A Padé approximant is the "best" approximation of the given function by a rational function of given order. For example, log 10 (95) can be accurately approximated by its Padé approximant at the expansion point x = 1, keeping in mind that log 10 (100) = 2 and using the fact that log 10 (95) = log 10 100 100 95 = log 10 (100) − log 10 100 100 95 ≈ 1.0526 is close to 1. The Padé approximant will be very accurate in this case, because 1.0526 is very close to the expansion point x = 1.
Evaluation of the decadic logarithm can be expressed by the natural logarithm log 10 (z) = ln(z) ln (10) , where ln(10) ≈ 2.302585093 is constant. The first task for students would be to make research related to Padé approximants. For students with interest in programming, the additional task would be to develop a code in appropriate software, such as Matlab. Padé approximants of different orders can be used for approximation of ln(z), for arguments z close to 1, i.e., z ≈ 1. As the expansion point z 0 = 1 is a root of ln(z), the accuracy of the Padé approximant decreases. Setting the "OrderMode" option in the Matlab Padé command to relative compensates for the decreased accuracy. Thus, here, the Padé approximant of the (m,n) order uses the form ln(z) ≈ , where α and β are coefficients (the coefficients of the polynomials need not be rational numbers). The Horner algorithm transforms polynomials into a computationally efficient form [30]. One possible procedure for solving the Colebrook equation using Padé approximants and the simple fixed-point iterative method is shown in Figure 2. Detailed step-by-step calculation using a table can be found in [30]. , where (10) ≈ 2.302585093 is constant. The first task for students would be to make research related to Padé approximants. For students with interest in programming, the additional task would be to develop a code in appropriate software, such as Matlab. Padé approximants of different orders can be used for approximation of ( ), for arguments z close to 1, i.e., ≈ 1. As the expansion point = 1 is a root of ( ), the accuracy of the Padé approximant decreases. Setting the "OrderMode" option in the Matlab Padé command to relative compensates for the decreased accuracy. Thus, here, the Padé approximant of the (m,n) order uses the form ( ) ≈ where α and β are coefficients (the coefficients of the polynomials need not be rational numbers). The Horner algorithm transforms polynomials into a computationally efficient form [30]. One possible procedure for solving the Colebrook equation using Padé approximants and the simple fixed-point iterative method is shown in Figure 2. Detailed step-by-step calculation using a table can be found in [30].  Students interested in programming can use tic and toc functions of Matlab and GNU/Octave [49,50] to measure the speed of calculations using the classical fixed-point method vs. the fixed-point method combined with Padé approximants (let us remind readers that GNU/Octave software is available for free [49]). Students can combine the Newton-Raphson method with Padé approximants and check the accuracy of Padé approximants of different orders [30].
For a general understanding of the concept of Padé approximants, 90 min is required (45 min for an introduction and a further 45 min for numerical tests). Then, that additional 90 min would be required for the application of Padé polynomials for solving the Colebrook equation (again, 45 min for an introduction and 45 min for numerical examples). Students with an interest in programming should have an additional 45-min exercise (a presentation of complexity vs. accuracy of Padé approximants).

Approximations by Multi-Point Methods with Internal Cycles
Some iterative methods belong to the group of multi-point methods [8] that are very powerful and that can reach the accurate solution of the Colebrook equation in the first iteration. An example of an approximation of the Colebrook equation based on internal cycles is given by Serghides [51]. Serghides' methods belong to Steffensen types of methods, which do not require derivatives (other Steffensen types of methods can be used in addition [52,53]).
In . Students should calculate this value and evaluate the introduced error compared with the accurate iterative solution [20]. An additional task could be to find additional approximations with internal cycles, such as Romeo et al. [55] or Offor and Alabi [56] and to discuss what additional strategies were used in these cases. For this task, 45 min is required.

Special Functions
The Colebrook equation is implicit with respect to the unknown flow friction factor. However, the Colebrook equation can be rearranged in the explicit form using the Lambert W function, where further, this special function can be evaluated only approximately [5]. Some numerical constraints in using the Lambert W function, which can occur in our case, will be detected [37] and a way to mitigate this inconvenience will be shown [27][28][29].
For tasks related to special functions, students need 45 min for a general introduction and an additional 90 min to work on examples.

Lambert W Function
Interested readers may find more information on the Lambert W function through [5,6,57,58]. Then, students should understand why some functions, such as logarithmic, square root, or exponential, are available in pocket calculators and spreadsheet solvers, while the Lambert W function is not [6].
In this task, students should understand how computers use a series expansion to interpret functions. Moreover, they need to understand how to use Visual Basic for Applications (VBA) modules in order to introduce new functions in MS Excel [59] (related to this task, students interested in programming should take at least 90 min).
The Lambert W function is used in mathematics to solve equations in which the unknown appears both outside and inside an exponential function or a logarithm, such as 3x + 2 = e x or x = ln(4·x). According to Fukushima [60], the Lambert W function is defined as the solution of a transcendental equation with κ = W(κ)·e W(κ) . Such equations cannot be, except in special cases, solved explicitly in terms of the finite number of arithmetic operations, nor in terms of exponentials or logarithms. In the case of the Colebrook equation, students can consult scientific literature, while here we will show formulation from [28,29], Equation (7), while a similar formulation can be found in [27]: Our carefully chosen examples, Re = 2.3 × 10 5 and ε D = 10 −4 , (2) Re = 4.6 × 10 7 and ε D = 0.037, respectively, are chosen to show that using floating-point arithmetic, Equation (7), gives the accurate solution for Example 1, while for Example 2, an arithmetic overflow error will occur. It is because W(e x ) becomes too large for ordinary computer registers.

Wright ω Function and Related Approximations
As shown in Biberg [27] and Brkić and Praks [28,29], the problem with too large numbers in computer registers can be solved by introducing a cognate of the Lambert W function-the Wright ω function. The Wright ω function y = ω(x) = W(e x ) is a solution to the equation y + ln(y) = x. Here, students should reformulate Equation (7) through the Wright ω function [50].
In Matlab, but also in the GNU/Octave software [49], students can use the free WrightOmegaq library [50]. In this case, the command is y = wrightOmegaq(210441.7).
Students can easily verify that the solution y satisfies the relation y + ln(y) = x. Let us reiterate that for Example 2, the value of x is approximately 210441.7.
In the following task, students should explore how to estimate the Wright ω function without Matlab or GNU/Octave. The simplification W(e x ) − x ≈ ln(x)· 1 x − 1 gives an accurate approximation suitable for engineering practice-Equation (8) [28,29]: where 2·2.51 ln(10) ≈ 2.18 and 2.18·3.71 ≈ 8.0878. Students should evaluate the relative error of the approximation, Equation (8), and develop a more accurate approximation using the enhanced model W(e x ) − x ≈ c + ln(x)· 1 x − 1 , where c is a constant. Of course, if c = 0, the enhanced model is equivalent to the model (8). For the general case, the constant c can be found by the MS Excel Solver, in order to minimize the relative error of the model with the iterative solution of the Colebrook equation provided by MS Excel. Students should do this optimization task with randomly selected input values of the Colebrook model Re, ε D , which are sampled from the domain of applicability of the Colebrook equation. Finally, students can compare the relative error of this enhanced model with the relative error of the original model.

Tania Function and Related Approximations
Students should examine the Tania function, which is defined as Tania(x) = ω(x + 1), Tania(x) = W(e x ), all with reference to [61] (the Euler T function should be an appropriate task to be further examined [39]).
For this task, 45-90 min should be assigned, during which students need to try to develop ability accepting new functions with new concepts [62].

Conclusions
Here we present how the relatively simple Colebrook equation, widely known by all students of hydraulics, can be used for introducing in curricula various iterative methods, special functions such as Lambert W, polynomials such as Padé approximants, optimization methods, etc.
The proposed exercises are suitable for students of hydraulics (civil, mechanical, chemical, or petroleum engineering curriculum, including the exploitation of oil and gas), but can be used as practical examples in teaching of numerical methods for students of mathematics. The proposed exercises can be extended for the teaching of pipe network analysis [31][32][33][34][35][36]. These methods can be combined with recent experimental works [63]. Furthermore, instead of pipes, a modified version of the Colebrook equation can be used for fuel cells [64,65].
Author Contributions: The paper is a product of the joint efforts of the authors, who worked together on models of natural gas distribution networks. P.P. has a scientific background in applied mathematics and programming, while D.B. has a background in control and applied computing in mechanical and petroleum engineering. They used the experience gained from their previous papers related to the Colebrook equation to develop proposed teaching exercises. D.B. wrote a draft of this article and prepared figures.