Evolutionary Optimization of Colebrook’s Turbulent Flow Friction Approximations

: This paper presents evolutionary optimization of explicit approximations of the empirical Colebrook’s equation that is used for the calculation of the turbulent friction factor ( λ ), i.e., for the calculation of turbulent hydraulic resistance in hydraulically smooth and rough pipes including the transient zone between them. The empirical Colebrook’s equation relates the unknown ﬂow friction factor ( λ ) with the known Reynolds number (R) and the known relative roughness of the inner pipe surface ( ε /D). It is implicit in the unknown friction factor ( λ ). The implicit Colebrook’s equation cannot be rearranged to derive the friction factor ( λ ) directly, and therefore, it can be solved only iteratively [ λ = f( λ , R, ε /D)] or using its explicit approximations [ λ ≈ f(R, ε /D)], which introduce certain error compared with the iterative solution. The optimization of explicit approximations of Colebrook’s equation is performed with the aim to improve their accuracy, and the proposed optimization strategy is demonstrated on a large number of explicit approximations published up to date where numerical values of the parameters in various existing approximations are changed (optimized) using genetic algorithms to reduce maximal relative error. After that improvement, the computational burden stays unchanged while the accuracy of approximations increases in some of the cases very signiﬁcantly.

The Colebrook equation is empirical, and hence, its accuracy can be disputed; it is still accepted in engineering practice as sufficiently accurate.It is still widely used in petroleum, mining, mechanical, civil and in all branches of engineering that deal with fluid flow.
Hydraulic resistance: Hydraulic resistance in general depends on the flow rate [53][54][55][56][57][58][59].To make things even more complex, hydraulic resistance is usually expressed through the flow friction factor, such as Darcy's (λ), where further pressure drop and flow rate are correlated with the well-known formula by Darcy and Weisbach.In the non-linear Darcy-Weisbach law for pipe flow, Darcy's friction Fluids 2017, 2, 15 2 of 27 factor (λ) is variable and always depends on flow.This assumption stands also if Fanning's friction is in use, since its physical meaning is equal to Darcy's friction (λ) with the difference that Fanning used the hydraulic radius while Darcy used the diameter to define the friction factor.Darcy's friction factor, known also under the names of Moody or Darcy-Weisbach, is four-times greater than Fanning's friction factor.The Fanning friction factor is more commonly used by chemical engineers and those following the British convention.
Colebrook equation: To be more complex, as already said, the widely-used empirical and nonlinear Colebrook's Equation (1) for calculation of Darcy's friction factor (λ) is iterative, i.e., implicit in the fluid flow friction factor since the unknown friction factor appears on the both sides of the equation [λ 0 = f(λ 0 , R, ε/D)] [2].This unknown friction factor (λ) cannot be extracted to be on the left side of the equal sign analytically, i.e., with no use of some kind of mathematical simplifications.Better to say, it can be expressed explicitly only if approximate calculus takes place.
λ 0 denotes the high precision iterative solution of Colebrook's equation, which is treated here as accurate; R denotes the Reynolds number; while ε/D denotes the relative roughness of the inner pipe surfaces.All three mentioned values are dimensionless.The Colebrook equation is also known as the Colebrook-White equation or simply the CW equation [1,2].This equation is valuable for the determination of hydraulic resistances for the turbulent regime in smooth and rough pipes including the turbulent zone between them, but it is not valid for the laminar regime.It describes a monotonic change in the friction factor (λ) during the turbulent flow in commercial pipes from smooth to fully rough.Moody's and Rouse's charts [3,4] represent the plots of the Colebrook equation over a very wide range of the Reynolds number (R from 2320 to 10 8 ) and relative roughness values (ε/D from 0 to 0.05).Besides some of its shortcomings [54], today, Colebrook's equation is accepted as the informal standard of accuracy for the calculation of the hydraulic friction factor (λ).
Accuracy: As already noted, the Colebrook equation is empirical, and therefore, its accuracy can be disputed; the equal sign "=" in "λ 0 = f(λ 0 , R, ε/D)", i.e., in Equation (1), instead of the approximately equal sign "≈", can be treated as accurate only conditionally [48].In this paper, the iterative solution of Colebrook's equation (λ 0 ) after enough iterations is treated as accurate and is used for comparison as the standard of accuracy where the accuracy of the friction factor (λ) calculated using the shown approximations will be compared with it.
Lambert W-function: The Colebrook equation can be rearranged in explicit form only approximately [λ ≈ f(R, ε/D)], where the approach with the Lambert W-function can be treated as a partial exemption from this rule [6][7][8][60][61][62], but also, further evaluation of the Lambert W-function function is approximate.
Looped network of pipes: The use of the accurate explicit approximations should be prioritized over the use of the iterative solution in the calculation of looped networks of pipes, since in that way, double iterative procedures, one for the Colebrook equation and one for the solution of the whole looped system of pipes, can be avoided [63][64][65][66][67].
Goal of the study: The goal is to increase the accuracy of the already available explicit approximation of Colebrook's equation.This is accomplished using genetic algorithms.

Genetic Algorithm Optimization Technique
Methodology: Genetic algorithms are one of the evolutionary computational intelligence techniques [45,46], inspired by Darwin's theory of biological evolution.Genetic algorithms provide solutions using randomly-generated strings (chromosomes) for different types of problems, searching the most suitable among chromosomes that make the population in the potential space of solutions.Genetic optimization is an alternative to the traditional optimal search approaches, which make it difficult to find the global optimum for nonlinear and multimodal optimization problems.Thus, genetic algorithms have been successful for example in solving combinatorial problems, control applications of parameter identification and control structure design, as well as in many other areas [47][48][49][50][51][52].
The used optimization approach: Here, the approach with genetic algorithms is implemented to optimize the parameters of the available approximations of the Colebrook equation for the hydraulic friction factor determination in order to improve their accuracy, at the same time retaining the previous complexity and computational burden of approximations.Small letters in the equations throughout the paper correspond to the numerical values before, while capital letters to the numerical values after optimization through genetic algorithms, as is picturesquely presented in Figure 1.
The used optimization approach: Here, the approach with genetic algorithms is implemented to optimize the parameters of the available approximations of the Colebrook equation for the hydraulic friction factor determination in order to improve their accuracy, at the same time retaining the previous complexity and computational burden of approximations.Small letters in the equations throughout the paper correspond to the numerical values before, while capital letters to the numerical values after optimization through genetic algorithms, as is picturesquely presented in Figure 1.Genetic algorithms are very powerful tools for optimization.Samadianfard [47] uses genetic programming, a sort of genetic algorithm, to develop his own explicit approximations to the Colebrook equation.Furthermore, genetic algorithms can be used together with some other techniques of artificial intelligence, such as neural networks [50][51][52].
Real coded genetic algorithms are used in this paper.The real coded genetic algorithms use the optimization designed cost function that minimizes maximal relative error, δmax, as follow (2): In (2), δ denotes relative (percentage) error, λ0 denotes the high precision iterative solution of Colebrook's equation, which is treated as accurate here, λ denotes the hydraulic friction factor solution calculated by each approximation considered and n denotes number of pairs of λ0 and λ used for optimization (in our case, n = 90,000).The fitness function is evaluated in a large number of 90 thousand points uniformly distributed in domains of the Reynolds number (R) and the relative roughness (ε/D).These domains correspond to those from the Moody diagram [3]; i.e., 2320 < Re < 10 8 , 10 6 < ε/D < 0.05, where roughness ε usually for PVC and plastic pipes is 0.0015-0.007mm, copper, lead, brass, aluminum (new) 0.001-0.002mm and for steel commercial pipe 0.045-0.09mm.The subjects of genetic optimization are coefficients in approximations, i.e., numeric coefficients in each approximation are changed by genetic algorithms in order to minimize the fitness function (2).In that way, approximations are changed in order to match the accuracy of the iterative solution of Colebrook's equation as much as possible.Simultaneous optimization of all coefficients in each approximation is attempted, while the range of values of parameters in which optimal solutions are searched is always in the arbitrary neighborhood of the initial values.Here, we chose to present the results obtained with the fitness function (2) in order to reduce the maximal error of each approximation as much as possible (assuming that the reduction of the maximal relative error is of the highest importance for the practical use of approximations).Genetic algorithms' performance depends on their parameter values, so genetic algorithm parameters were carefully selected by conducting numerous experiments.In the implemented algorithm, a real-coded population of 100 individuals, an elitism of 10 individuals and a scattered crossover function are used.All of the members are subjected to adaptive feasible mutation except for the elite.The individuals are randomly selected by the Genetic algorithms are very powerful tools for optimization.Samadianfard [47] uses genetic programming, a sort of genetic algorithm, to develop his own explicit approximations to the Colebrook equation.Furthermore, genetic algorithms can be used together with some other techniques of artificial intelligence, such as neural networks [50][51][52].
Real coded genetic algorithms are used in this paper.The real coded genetic algorithms use the optimization designed cost function that minimizes maximal relative error, δ max , as follow (2): In (2), δ denotes relative (percentage) error, λ 0 denotes the high precision iterative solution of Colebrook's equation, which is treated as accurate here, λ denotes the hydraulic friction factor solution calculated by each approximation considered and n denotes number of pairs of λ 0 and λ used for optimization (in our case, n = 90,000).The fitness function is evaluated in a large number of 90 thousand points uniformly distributed in domains of the Reynolds number (R) and the relative roughness (ε/D).These domains correspond to those from the Moody diagram [3]; i.e., 2320 < Re < 10 8 , 10 6 < ε/D < 0.05, where roughness ε usually for PVC and plastic pipes is 0.0015-0.007mm, copper, lead, brass, aluminum (new) 0.001-0.002mm and for steel commercial pipe 0.045-0.09mm.The subjects of genetic optimization are coefficients in approximations, i.e., numeric coefficients in each approximation are changed by genetic algorithms in order to minimize the fitness function (2).In that way, approximations are changed in order to match the accuracy of the iterative solution of Colebrook's equation as much as possible.Simultaneous optimization of all coefficients in each approximation is attempted, while the range of values of parameters in which optimal solutions are searched is always in the arbitrary neighborhood of the initial values.Here, we chose to present the results obtained with the fitness function (2) in order to reduce the maximal error of each approximation as much as possible (assuming that the reduction of the maximal relative error is of the highest importance for the practical use of approximations).Genetic algorithms' performance depends on their parameter values, so genetic algorithm parameters were carefully selected by conducting numerous experiments.In the implemented algorithm, a real-coded population of 100 individuals, an elitism of 10 individuals and a scattered crossover function are used.All of the members are subjected to adaptive feasible mutation except for the elite.The individuals are randomly selected by the Roulette method.Optimization with genetic algorithms was carried out in MATLAB R2010a by MathWorks (Natick, MA, USA).The practical domain of the Reynolds number (R) and relative roughness of inner pipe surface (ε/D) are covered by a mesh of n = 90,000 points for this optimization.In these 90 thousand points, the iterative solution of the implicitly-given Colebrook equation, λ 0 , and the non-iterative solution for every single observed approximation, λ, are calculated.The optimization of every single approximation lasts several hours.All evaluations of error were performed in MATLAB, with further confirmations in MS Excel to maintain full comparability with the study of Brkić [10] (for the use of iterative calculus in MS Excel Ver.2007, see Brkić [11]; in Brkić and Tanasković [68], MS Excel is also used for other extensive, but non-iterative calculations).The mesh in MS Excel over the practical domain of the Reynolds number (R) and relative roughness of the inner pipe surface (ε/D) consists of n = 740 uniformly-distributed points.
Alternative optimization approaches: The main goal of the optimization in our case is to reduce the maximal error (δ max ) of the every single observed approximation.This means that sometimes, the average (mean) relative error in the practical range of the Reynolds number (R) and the relative roughness of inner pipe surface (ε/D) increases compared to the model of the observed approximation with initial, non-optimized values of the parameters.Of course, using genetic algorithm optimization with the function defined to reduce maximal error, this error is reduced more or less efficiently, which at the same time does not mean that average error is necessarily increased or decreased.Although the minimization of average error is not set as a goal by Equation ( 2), it can be reduced also during the optimization.Instead of the here already shown fitness function Equation (2), it can be redefined to simultaneously reduce average and maximal error Equation (3).In that way, both errors, i.e., maximal relative error and average (mean) relative error, can be reduced simultaneously for sure.This requires more one-off computational efforts compared with the approach in which only one type of error is reduced; in our case, maximal relative error, δ max , while the fitness function is defined by Equation (2).In Equation (3), the first term reduces average (mean) relative error, δ avr ; the second term reduces maximal error δ; while weights k 1 and k 2 can be used to signify one of the terms and reduce the influence of other.In that case, a compromise between the reduction of the maximal and average relative error is obtained.
Furthermore, the fitness function can be set to reduce simultaneously the mean square error δ MSE and maximal relative error δ, as in Equation (4).As already noted for Equation (3), the ratio between weight coefficients k 3 and k 4 determines the influence of the mean square error δ MSE and the maximal relative error δ in the optimization.According to many different criteria, the values of coefficients in the existing explicit approximation to the Colebrook equation can be used.Using many computational resources, all three errors shown in our paper can be simultaneously reduced Equation ( 5), but such a procedure seems to be quite elusive.

Explicit Approximations of Colebrook's Equation
Colebrook's Equation (1) [2] suffers from being implicit in the unknown friction factor (λ).It requires an iterative solution where convergence to the final accuracy of the observed approximation typically requires less than seven iterations.As Brkić [10] proposed, here is used even a few thousand iterations to be sure that a sufficient value of accuracy for the friction factor, λ 0 , is reached.
As already stated, the implicit Colebrook's equation cannot be rearranged to derive the friction factor directly in one step, while iterative calculus can cause a problem in the simulation of flow in a pipe system in which it may be necessary to evaluate the friction factor hundreds or thousands of times.This is the main reason for attempting to develop a relationship that is a reasonable and an as accurate as possible approximation for the Colebrook equation, but which is explicit in the friction factor.These approximations are used for the calculation of the friction factor (λ), which is compared with the very accurate solution (λ 0 ) calculated using the iterative procedure.
The accuracy of existing approximations of Colebrook's equation was thoroughly checked by many researchers [10][11][12][13][14][15][16].Yıldırım [14] conducted a comprehensive analysis of existing correlations for single-phase friction, but he used Techdig 2.0 software to read the date from the Moody diagram, which caused remarkable reading error.One must be always aware that the Moody diagram [3] was constructed using Colebrook's equation [2] and not opposite.After all, the main conclusion of all papers [10][11][12][13][14][15][16] is that the relative error, δ, is non-uniformly distributed over the domain of the Reynolds number (R) and the relative roughness (ε/D).
The relative error δ is defined in Equations ( 2)-( 4) of this paper, the average (mean) relative error δ avr in Equation ( 3) and the mean square error δ MSE in Equation (4).All three types of error are used in further text for the estimation of the accuracy of the examined explicit approximations of the Colebrook equation, but the accent is on the minimization of the maximal relative error, δ max .
Using the shown genetic algorithm optimization technique, the values of existing parameters of the explicit approximations are improved compared to the iterative solution of Colebrook's equation.This means that the error of approximations decreases while the computational burden stays unchanged.In this section, new parameters are shown, and the reduction of the maximal relative error, δ max , is estimated.The relative errors of the approximations shown in further text of this paper are calculated as δ = [(λ − λ 0 )/λ 0 ]•100%, where λ is the Darcy friction factor calculated using the observed approximation, while λ 0 is the iterative solution of Colebrook's equation, which can be used as accurate after enough iterations (here, set to the maximal available number of iterations in MS Excel, which is 32,767, as explained in [10]).
Each of the 25 observed approximations is supplied with three diagrams; the first is the distribution of the relative error over the practical domain of applicability in engineering practice; the second is the same as the first, but with the relative error distribution after optimization; and the third is a comparative diagram.For the first two mentioned diagrams, the entire practical domain of the Reynolds number (R) and the relative roughness of the inner pipe surface (ε/D) is covered with a 740 point-mesh (diagrams produced in MS Excel).For the first two figures with approximations, same Fluids 2017, 2, 15 6 of 27 pace of error is used for non-optimized and for the optimized approximation, to provide a more easy comparison with the exceptions of Approximation (Appr.)10; Equation (15), Appr.11; Equation ( 16) and Appr.14; Equation (19), where the optimization was extremely successfully performed.The mesh of 740 points is formed in MS Excel using 20 values of the relative roughness (ε/D) (shown in the related figures) and using 37 values of the Reynolds number (R); from 10 4 -10 5 with a pace of 10 4 , from 10 5 -10 6 with a pace 10 5 , from 10 6 -10 7 with a pace of 10 6 and from 10 7 -10 8 with a pace of 10 7 .
According to Winning and Coole [16], using the value of mean square error δ MSE defined in ( 4), all approximations can be classified into four groups (the very small error is lower than 10 −11 ; small is between 10 −11 and 10 −8 ; medium is between 10 −8 and 5 × 10 −6 ; and large is above 5 × 10 −6 ).This criterion is also used in further evaluation.
Regarding accuracy, it should be noted that the inner roughness of the pipe, ε, cannot be determined easily [17], so the physical interpretation of the relative roughness of the inner pipe surface (ε/D) is not the subject of this study.
For genetic algorithm optimization, MATLAB 2010a by MathWorks is used.For this purpose, a mesh of 90 thousand points over the entire practical domain of the Reynolds number (R) and the relative roughness of the inner pipe surface (ε/D) is generated.For these 90 thousand pairs of Reynolds number (R) and the relative roughness of inner pipe surface (ε/D), the friction factor (λ 0 ) is very accurately calculated to be used as a pattern during the procedure of optimization.Although genetic optimization allows for virtually all coefficients in each and every formula considered to be optimized in the search for the best result, which can be considered as an important advantage of the presented technique, the selection of the coefficients included in the optimization for each formula was the object of careful consideration.Not only the inclusion of more coefficients significantly widens the search space to be covered, it appears also that some approximations are highly sensitive to the changes of some coefficients.Therefore, subsets of coefficients to be optimized within each formula were selected on the basis of multiple trials, approaches considered by other authors and our own vast experience.The number of digits in each optimized coefficient was also limited as a tradeoff of the desired punctuality of approximation and the practical usability of the formula.
The efficiency of computing in the computer environment stays unchanged between non-optimized and related optimized approximations, since the model of the approximation stays unchanged; i.e., the number of logarithmic and power expressions stays unchanged [9,18].Only the change of integer power to non-integer power in some approximation can increase the computational burden, but even then not significantly.
In the following Figures 2-26, symbols and zones with green and red color represent: the ∆δ-decreased level of maximal relative error δ max ; (1) zone of increased relative error δ (red); (2) zone of decreased relative error δ (green).

Conclusions
Today, Colebrook's equation is mostly accepted as an informal standard for the modeling of turbulent flow in hydraulically smooth and rough pipes, including the transient zone in between.Of course, approximations carry certain error compared with the iterative solution where the highest level of accuracy can be reached after enough iterations.The explicit approximations give a relatively good prediction of the friction factor (λ) and can reproduce accurately Colebrook's equation and its Moody's plot.Usually, more complex models of approximations are more accurate and vice versa.Using genetic algorithms in order to increase the accuracy of available approximations of the Colebrook equation for flow friction, the numerical values of empirical parameters in 25 existing models of approximations are changed while the computational burden remains the same.Using the value of decreased maximal relative error, Δδ, and the change of relative error over the entire domain of the Reynolds number (R) and the relative roughness of inner pipe surface (ε/D), the success of genetic optimization is summarized in Table 1.Brkić approximation (Appr.2): Relevant parameters and errors related to the approximation by Brkić [19] after and before optimization (7) (Appr.2) are given in Figure 3.
Sonnad and Goudar approximation (Appr.10): Relevant parameters and errors related to the approximation by Sonnad and Goudar [26] after and before optimization (15) (Appr.10) are given in Figure 11.Vatankhah and Kouchakzadeh [43,44]    Sonnad and Goudar approximation (Appr.10): Relevant parameters and errors related to the approximation by Sonnad and Goudar [26] after and before optimization (15) (Appr.10) are given in Figure 11.Vatankhah and Kouchakzadeh [43,44] by changing parameter ɑ11 to A11-0.31 slightly change the model of Sonnad and Goudar [26].They used line fitting tool for optimization.We failed with further optimization using genetic algorithms. ( 11 to A 11 -0.31 slightly change the model of Sonnad and Goudar [26].They used line fitting tool for optimization.We failed with further optimization using genetic algorithms. Romeo, Royo and Monzón approximation (Appr.11): Relevant parameters and errors related to the approximation by Romeo et al. [27] after and before optimization (16) (Appr.11) are given in Figure 12; already shown in the form of the preliminary note in Ćojbašić and Brkić [42].
Chen approximation (Appr.13): Relevant parameters and errors related to the approximation by Chen [29] after and before optimization (18) (Appr.13) are given in Figure 14.
Serghides approximation (Appr.14): Relevant parameters and errors related to the approximation by Serghides [30] after and before optimization (19) (Appr.14) are given in Figure 15.This optimization is already shown in the form of the preliminary note in Ćojbašić and Brkić [42].
Zigrang and Sylvester approximation (Appr.17): Relevant parameters and errors related to the approximation by Zigrang and Sylvester [32] after and before optimization (22) (Appr.17) are given in Figure 18.
Barr approximation (Appr.19): Relevant parameters and errors related to the approximation by Barr [33] after and before optimization (24) (Appr.19) are given in Figure 20.
Round approximation (Appr.20): Relevant parameters and errors related to the approximation by Round [34] after and before optimization (25) (Appr.20) are given in Figure 21.Chen approximation (Appr.21) Relevant parameters and errors related to the approximation by Chen [36] after and before optimization (26) (Appr.21) are given in Figure 22.Eck approximation (Appr.23): Relevant parameters and errors related to the approximation by Eck [38] after and before optimization (28) (Appr.23) are given in Figure 24.
Wood approximation (Appr.24): Relevant parameters and errors related to the approximation by Wood [39] after and before optimization (29) (Appr.24) are given in Figure 25.
Moody approximation (Appr.25): Relevant parameters and errors related to the approximation by Moody [40] after and before optimization (30) (Appr.25) are given in Figure 26.

Conclusions
Today, Colebrook's equation is mostly accepted as an informal standard for the modeling of turbulent flow in hydraulically smooth and rough pipes, including the transient zone in between.Of course, approximations carry certain error compared with the iterative solution where the highest level of accuracy can be reached after enough iterations.The explicit approximations give a relatively good prediction of the friction factor (λ) and can reproduce accurately Colebrook's equation and its Moody's plot.Usually, more complex models of approximations are more accurate and vice versa.Using genetic algorithms in order to increase the accuracy of available approximations of the Colebrook equation for flow friction, the numerical values of empirical parameters in 25 existing models of approximations are changed while the computational burden remains the same.Using the value of decreased maximal relative error, ∆δ, and the change of relative error over the entire domain of the Reynolds number (R) and the relative roughness of inner pipe surface (ε/D), the success of genetic optimization is summarized in Table 1.Since the main idea of this study was to use metaheuristic optimization to improve the accuracy of a wider set of Colebrook's turbulent flow friction approximations, as opposed to other approaches used by other authors and ourselves, Genetic Algorithms were selected as being referential among global optimization techniques.Additional accuracy for the certain approximations can be possibly reached through rearrangement of their structure (such as was done for Appr. 10 [26,43]), using the MS Excel fitting tool [49], etc.Furthermore, some further use of genetic algorithms can be encouraged.Future
Swamee and Jain approximation (Appr.22): Relevant parameters and errors related to the approximation by Swamee and Jain [37] after and before optimization (27) (Appr.22) are given in Figure 23.
Swamee and Jain approximation (Appr.22): Relevant parameters and errors related to the approximation by Swamee and Jain [37] after and before optimization (27) (Appr.22) are given in Figure 23.

Table 1 .
Maximal relative error of the explicit approximations of the Colebrook-White equation before and after genetic optimization.
by changing parameter

Table 1 .
Maximal relative error of the explicit approximations of the Colebrook-White equation before and after genetic optimization.