Accurate and Efficient Explicit Approximations of the Colebrook Flow Friction Equation Based on the Wright ω-Function: Reply to the Discussion by Majid Niazkar

Abstract

In this reply, we present updated approximations to the Colebrook equation for flow friction. The equations are equally computational simple, but with increased accuracy thanks to the optimization procedure, which was proposed by the discusser, Dr. Majid Niazkar. Our large-scale quasi-Monte Carlo verifications confirm that the here presented novel optimized numerical parameters further significantly increase accuracy of the estimated flow friction factor.

1. Introduction

We thank discusser Dr. Majid Niazkar for his valuable discussion of our article [1]. In his discussion, he proposed to optimize numerical parameters of our explicit approximations of the Colebrook’s equation, in order to increase their accuracy in comparison to the implicitly given Colebrook’s flow friction factor [2]. The discusser also presented optimized numerical parameters for the appropriate segment of the Wright ω-function which is used in our approximations of the Colebrook equation.

The original Colebrook equation for flow friction [3] is of implicit nature, and hence our goal is to develop its computational simple but accurate explicit approximations. It is also of empirical nature [4] and hence its accuracy can be disputed [5], but equally because of engineering comparisons, the developed approximations should be as accurate as possible.

Niazkar uses ten criterions to evaluate accuracy of the proposed equations (six local evaluation criterions); Niazkar’s Equations (4)–(9) and four global; Niazkar’s Equations (10)–(13), but here we are focused on minimization of the maximal absolute relative error in percentages (basically referred to Niazkar’s local criterions given by his Equations (8) and (9)). The equations proposed by Niazkar have the same structure as ours, i.e., they follow our models, but their accuracy increases by changing numerical values through optimization procedures (Niazkar’s Model 3 improves our Model 1, while his Model 4 improves our Model 2; Models 1 and 2 are Equations (5) and (6) of [1], respectively). In that way, computational burden related to these models remains unchanged while accuracy increases.

Accuracy of our equations from [1] (Models 1 and 2 in Niazkar’s Discussion [2]) has been already confirmed in [6,7]. Using the methodology proposed by Niazkar [8], our approximation can be classified in the group of one-step relations which clearly confirm their simplicity.

In this reply we present updated numerical parameters for these models, which are equally computational simple, but with further increased accuracy. We highly appreciate Niazkar’s approach to optimize numerical parameters of our approximations and we can confirm its efficiency in error reduction. Comparable approaches can be found in the scientific literature [9,10,11], such as optimization of the other already available approximations of the Colebrook equations using genetics algorithms, etc.

2. Updated Approximations

Approximations given in Equation (3) of our discussed article [1] are based on the approximation of the Wright ω-function as given in [12], while Equations (5) and (6) of [1] are found using symbolic regression software Eureqa [13,14] and they are also based on the approximation of the same branch of the Wright ω-function as in Equation 3. Here we will show approximations with the same structure as in [1], but with optimized numerical values of the parameters of the appropriate equations which give reduced absolute maximal relative errors. This optimization procedure significantly reduces the maximal relative error of approximations from our original paper. For Equation (3) of [1], the maximal relative error was reduced by 50%, and for Equation (5) of [1] by 42%, while for Equation (6) of [1] is improved only by 3% as the original model was close to the optimum. We used the local “fminsearch” nonlinear solver, and then the global optimization procedure “globalsearch”, both built-in Matlab optimization tools [15]. We reduced also the number of required digits for optimized constants by numerical experiments, in order to provide simple, but still accurate flow friction approximations. All the here presented approximations over-perform our original article [1], thanks to the global optimization of parameters [16].

In all our experiments, eight million quasi Monte-Carlo pairs [17] cover thoroughly the domain of the Reynolds numbers $Re$ and the relative roughness $\epsilon$ which is interest for engineering practice (all details can be found in [1]) and all numerical values of our approximations are optimized in that domain; 4000 < $Re$ < 10⁸ and 0 < $\epsilon$ < 0.05. In addition, not only empirical numerical values, but also those of clearly defined parameters, such as 0.8686 which is obtained from $\frac{2}{\ln \left(10\right)}$, are also optimized.

The general exact explicit formula of the solution of the Colebrook equation using the Wright ω-function with the vector of parameters $p=\left(\frac{2}{ln\left(10\right)}, ln\left(\frac{2·2.51}{\ln \left(10\right)}\right),\frac{2·2.51·3.71}{ln\left(10\right)}\right)$ can be implemented in Matlab through the “wrightOmegaq” Function [18] as given with the following code:

$$\begin{array}{c}B=log\left(R\right)-p\left(2\right);\\ A=R.\ast K./p\left(3\right);\\ \begin{array}{c}x=A+B;\\ f=1./\left(p\left(1\right)\ast \left(B+y\right)\right).^2\end{array}\end{array}{\displaystyle \rangle }$$

where $y=wrightOmegaq\left(x\right)-x$, $R$ is the vector of the dimensionless Reynolds number; $Re$ in Equation (1), $K$ is the vector of the dimensionless relative roughness of inner pipe surface; $\epsilon$ in Equation (1), $f$ is the dimensionless Darcy–Weisbach flow friction factor; $f$ in Equation (1), while the original vector of the approximate parameters [1] used for optimization is $p=\left(0.8686, 0.7794, 8.0878\right)$. In Matlab, “$log$” denotes natural logarithm. This Matlab code corresponds to Equation (1):

$$\begin{array}{c}f=\frac{1}{{\left(0.8686·\left(B+y\right)\right)}^2}\\ A\approx \frac{Re·\epsilon }{8.0878}\\ \begin{array}{c}B\approx \ln \left(Re\right)-0.7794\\ x=A+B\\ y=\omega \left(x\right)-x\end{array}\end{array}\},$$

(1)

2.1. Optimized Equation (3) of [1]

For optimization of Equation (3) of [1], we again used “fminsearch” and then “globalsearch” optimization procedures of Matlab. Approximation of the y has form $y \approx \frac{\ln \left(x\right)}{x}-\ln \left(x\right)$, while $y=wrightOmegaq\left(x\right)-x$ is treated as accurate. The optimized vector is now $p=\left(0.86902384, 0.7829415, 8.11718121\right)$ which reduces the maximal absolute relative error of Equation (3) of [1] from 0.1523% to 0.100793%– 0.101% for the updated Equation (3) of [1] which is here given in Equation (2) (from the plot in Figure 1, we can see how the maximal relative error was reducing during the optimization process). The number of relevant digits in Equation (2) was found by numerical experiments.

$$\begin{array}{c}f=\frac{1}{{\left(0.86902384·\left(B+y\right)\right)}^2}\\ A\approx \frac{Re·\epsilon }{8.11718121}\\ \begin{array}{c}B\approx \ln \left(Re\right)-0.7829415\\ x=A+B\\ y\approx \ln \left(x\right)·\left(\frac{1}{x}-1\right)\end{array}\end{array}\},$$

(2)

Thus, this procedure reduces the maximal absolute relative error of Equation (3) of [1] from 0.1523% to 0.100793% while computation burden remains unchanged.

2.2. Optimized Equation (5) of [1]

The original Equation (5) of [1] in general can be rewritten as Equation (3).

$$\begin{array}{c}f=\frac{1}{{\left(p_1·\left(B+y\right)\right)}^2}\\ A=\frac{Re·\epsilon }{p_3}\\ \begin{array}{c}B=\ln \left(Re\right)-p_3\\ x=A+B\\ y\approx \ln \left(x\right)·\left(\frac{p_4}{x+p_5}-1\right)\end{array}\end{array}\},$$

(3)

where the original parameters $(p_1, p_2,p_3, p_4,p_5)$ in Equation (5) of [1] are given by the vector $p=\left(0.8686, 0.7794, 8.0878, 1.038, 0.332\right)$, while Niazkar’s proposal for the updated vector is $p=\left(0.86855,\ln \left(2.18\right), 8.0878, 1.03891, 0.33623\right)$. Niazkar refers in his Discussion to Equation (5) of [1], i.e., to Equation (3) with the original vector $p$ as Model 1, while with the vector $p$ updated by him as Model 3 [2]. The structure of the approximation of $\omega \left(x\right)-x$ for the observed case remains $\frac{p_4·\ln \left(x\right)}{x+p_5}-\ln \left(x\right)$. In Niazkar’s discussion, the symbol B is defined as $B\approx \frac{\ln \left(Re\right)}{2.18}$ while we originally proposed simpler version as $B\approx \ln \left(Re\right)-0.7794$ which is faster, as it does not require a division. The original version from [1], i.e., Model 1 as referred in Niazkar’s Discussion [2] has maximal relative error of no more than 0.0522%, while Niazkar’s Model 3 reduces the maximal relative error of which is now up to 0.0459%. As there is only a change of constants of the model, the complexity of the approximation is the same: our experiments show that speed of the optimized approximation remains unchanged (0.41 s for 8 million of Quasi Monte-Carlo pairs).

Further, for the purpose of this reply, we will optimize constants of Equation (5) of [1]. In other words, we will update all five constants of the model given with vector $p$ of Equation (3).

Consequently, the explicit approximation of the Colebrook following equation Equation (5) of [1] with the vector of parameter $p$, can be written in Matlab by the following code:

$$\begin{array}{c}B=log\left(R\right)-p\left(2\right);\\ A=R.\ast K./p\left(3\right);\\ \begin{array}{c}x=A+B;\\ f=1./\left(p\left(1\right)\ast \left(B+y\right)\right).^2;\\ y=\left(p\left(4\right) \ast log\left(x\right)\right)./\left(x+p\left(5\right)\right)-log\left(x\right)\end{array}\end{array}{\displaystyle \rangle },$$

As we already mentioned, symbol log(x) represents in the Matlab’s code the natural logarithm.

The original approximation given by Equation (5) of [1] can be coded by the vector $p=\left(0.8686, 0.7794, 8.0878, 1.038, 0.332\right)$. In our optimization, we minimized the maximal relative error of the flow friction factor $f$ by changing the vector $p$ (in contrary, Niazkar analyzed many types of errors in his Discussion [2]) and finally we obtained the following optimal parameters $p=\left(0.868585, 0.78157, 8.099752, 1.04796, 0.36322\right)$. The absolute maximal relative error is now no larger than 0.0366%.

Thus, this procedure reduces the maximal absolute relative error of Equation (5) of [1] from 0.0522% to 0.0366% (while in Niazkar’s Discussion it is up to 0.0459% [2]).

2.3. Optimized Equation (6) of [1]

Equation (6) of [1] is Model 2 of Niazkar’s Discussion [2], while his related version is Model 4. This approximation is given here as Equation (4).

$$\begin{array}{c}f=\frac{1}{{\left(p_1·\left(B+y\right)\right)}^2}\\ A\approx \frac{Re·\epsilon }{p_3}\\ \begin{array}{c}B\approx \ln \left(Re\right)-p_2\\ x=A+B\\ \begin{array}{c}C=\ln \left(x\right)\\ y\approx \frac{p_4·C}{x}+\frac{C-p_5}{x^2}-C\end{array}\end{array}\end{array}\},$$

(4)

The vector with the original parameters from [1], $p=\left(0.8686, 0.7794, 8.0878, 1.0119, 2.3849\right)$ gives an approximation with the maximal relative error of no more than 0.00845%. Niazkar’s vector $p=\left(0.86859,\ln \left(2.18\right), 8.0878, 1.01151, 2.37718\right)$ gives an approximation with the slightly reduced error. Using the optimization approach as described in Section 2, the approximation of $\omega \left(x\right)-x$ can be expressed as $y\approx \frac{p_4·\ln \left(x\right)}{x}+\frac{\ln \left(x\right)-p_5}{x^2}-\ln \left(x\right)$. We get the optimized numerical values of the vector as $p=\left(0.868558, 0.77898, 8.0861744, 1.011746, 2.3872\right)$ which reduces the maximal relative error which is now up to 0.00807% (the original approximation with the original numerical values from [1] is almost optimal and so reduction of the error is minimal). In contrary to previous cases, the optimized procedure was less effective. Optimization process of “fminsearch” was monitored and the related plot can be seen in Figure 2.

Thus, this procedure reduces the maximal absolute relative error of Equation (5) of [1] from 0.00845% to 0.00807592%.

3. Conclusions

In this reply, we verified the optimization of numerical coefficients of the explicit approximations of the Colebrook equation based on various approximations of the Wright ω-function. Our large-scale quasi Monte-Carlo simulations confirm that the optimization approach suggested by the discusser Dr. Majid Niazkar can reduce the maximal relative error of approximations. We also presented three updated approximations with optimized coefficients.

The presented results are based on our previous works [1,9,19,20,21,22,23,24,25,26,27].