$$ f^{-1/2} = -2 \text{log}_{10} \left( \frac{Rr}{3.7} + \frac{2.51}{Re} f^{-1/2} \right) $$
If you plot both sides of this, you obtain the following.
Simple bisection method would work just fine to solve this, but there's one detail missing - what do you choose as the bounds. I'll call the lower and upper bounds $f_{min}$ and $f_{max}$. You could say "for problems I do the friction factor is usually between 0.0wy and 0.0xz", but that's somewhat inelegant, and it risks failing without notice.
Robust lower and upper bounds
My solution is to derive an expression for the lower and upper bound for the friction factor. I will do this directly from the mathematics. Refer back to the graph I presented. The specifics are unimportant. What's important is the fact that no matter what the values in the equation ($Re$, $Rr$, and so on), the two lines will still be sloping the same way. This is the way in for me to formulate an inequality for the lower bound and upper bound.
Lower bound
It is absolutely true that for any value of f that is less than the intersection, the left hand side (LHS) will be greater than the right hand side (RHS). Echoing the equation with a new variable for the lower bound, $f_{min}$, gives:
$$ f^{-1/2} > -2 \text{log}_{10} \left( \frac{Rr}{3.7} + \frac{2.51}{Re} f^{-1/2} \right) $$
Thus, any expression for $f_{min}$ that satisfies this inequality under all conditions is absolutely satisfactory as a minimum for doing bisection method. Since any value will do (as long as it meets the condition), we can use some tricks. Firstly, note that the LHS is positive. Since any $f$ or $f_{min}$ is positive, one divided by the square root will also be positive. That means that if we just... set that equal to zero, the inequality is still valid to find a left-bounding value for $f_{min}$.
$$ 0 > -2 \text{log}_{10} \left( \frac{Rr}{3.7} + \frac{2.51}{Re} f^{-1/2} \right) $$
Indeed, if the previous inequality is right, this inequality is more right. That has utility because now we can solve this. For the log to be positive, the argument has to be greater than 1.
$$ \frac{2.51}{Re} f_{min}^{-1/2} > 1 - \frac{Rr}{3.7} $$
$$ f_{min}^{-1/2} > \frac{Re}{2.51} \left( 1 - \frac{Rr}{3.7} \right) $$
$$ f_{min} > \left( \frac{2.51}{Re} \right)^2 \left( 1 - \frac{Rr}{3.7} \right)^{-2} $$
The inequality sign is a little confusing here, but what that's really saying is that the value of f must be greater than that expression. So the expression establishes our lower bound.
Upper bound
EDIT:
A subsequent post supersedes this one:
http://gravitationalballoon.blogspot.com/2013/03/a-robust-method-for-numerical-solution.html
The solution presented here works but that one is much better. I will still leave the work here because it goes into much greater depth.
-----------------------
Similarly, by graph inspection we determine that for an upper bound, $f_{max}$, we will seek to satisfy:
$$ f_{max}^{-1/2} < -2 \text{log}_{10} \left( \frac{Rr}{3.7} + \frac{2.51}{Re} f_{max}^{-1/2} \right) $$
I'll use some simple algebra. I also have a specific form in mind that I'm looking for. If I can isolate a $\log \left( 1+x \right)$ term, then provided the $x$ is positive, the term will be positive, and in these inequalities, that's a opportunistic point for simplification.
$$ f_{max}^{-1/2} < -2 \text{log}_{10} \frac{Rr}{3.7} -2 \text{log}_{10} \left( 1 + \frac{3.7}{Rr} \frac{2.51}{Re} f_{max}^{-1/2} \right) $$
This is difficult because we can't use the same approach as the lower bound. Observe that removing the LHS or the term on the far right in isolation won't be valid. Doing so would make the inequality no longer universally true. We're allowed to add more to the LHS or take away more from the RHS, but not the other way around. My solution is thus the following trick:
$$ \ln \left( 1 + x \right) < x $$
This is true for all values of x greater than zero (or any real value). Use this to substitute into the previous inequality in place of the RHS log term.
$$ f_{max}^{-1/2} < -2 \text{log}_{10} \frac{Rr}{3.7} - \frac{2}{\ln{10}} \frac{3.7}{Rr} \frac{2.51}{Re} f_{max}^{-1/2} $$
Now just do algebra to get to the final form.
$$ \left( 1 + \frac{2 \times 3.7 \times 2.51 }{Rr Re \ln{10}} \right) f_{max}^{-1/2} < -2 \text{log}_{10} \frac{Rr}{3.7} $$
The relative roughness really should be less than 1, and we should change the sign to reflect this, particularly dealing with inequalities.
$$ \left( 1 + \frac{2 \times 3.7 \times 2.51 }{Rr Re \ln{10}} \right) < 2 \text{log}_{10} \frac{3.7}{Rr} \sqrt{ f_{max} } $$
$$ \left( \frac{ 1 + \frac{2 \times 3.7 \times 2.51 }{Rr Re \ln{10}} }{2 \text{log}_{10} \frac{3.7}{Rr} } \right) < \sqrt{ f_{max} }$$
$$ f_{max} > \left( \frac{ 1 + \frac{2 \times 3.7 \times 2.51 }{Rr Re \ln{10}} }{2 \text{log}_{10} \frac{3.7}{Rr} } \right)^2 $$
Like the last case, what this is really saying is that the true value of f is less than the expression that's been obtained here.
Special case of zero roughness
The Colebrook-White formula is still perfectly valid if the roughness value is actually zero. Unfortunately, one look shows that this is incompatible with at least one of the formulas above. If you plug in Rr=0 in the original equation, then you can run through the calculations to get the same thing for fmin, but not fmax. In this special case of fmax with zero roughness, the inequality we seek to identify is.
$$ f_{max}^{-1/2} < -2 \text{log}_{10} \left( \frac{2.51}{Re} f_{max}^{-1/2} \right) $$
This produces a uniquely different problem. The same general approach still works. In the approach I used I made a different substution, utilizing the fact that -1/sqrt(x) < ln(x), which is a slightly different variation on what was used above. Carrying it through to its conclusion gives the following alternative expression.
$$ f_{max} > \left( \frac{ \ln{10} + 1 }{ 2 \ln{ \frac{Re}{2.51} } } \right)^2 $$
The similarity to the previous fmax is obvious, but here we've avoided any potential hazard of dividing by zero.
Code for VBA
To solve this with a simple bisection method, it's probably best to formulate it as a root-finding problem. So say that 0 = RHS - LHS. One small note of caution is that the function log() means natural log in Visual Basic. In formulas it's base 10 log. Is this blatantly inconsistent? Yes. I often make it my policy to only use natural log in order to avoid confusion, which is hard to do.
$$ g(f) = -2 \text{log}_{10} \left( \frac{Rr}{3.7} + \frac{2.51}{Re} f^{-1/2} \right) - f^{-1/2} $$
When you actually type out the function it looks like -2 / Log(10) * Log(Rr / 3.7 + 2.51 / Re * f ^ (-1 / 2)) - f ^ (-1 / 2). For argument ordering, I'm using Rr followed by Re, to be consistent with the usual way of writing the equation. Here is the code, which is in working order in my sheet:
Function Colebrook(Rr As Double, Re As Double) As DoubleTo use this, if you're in Excel 2007, you need to have your workbook saved as the xlsm file type. Then go to the Developer tab, and click the Visual Basic function. In the window on the left you need to right-click the VBAProject + other stuff name and Insert. Then you should be good to paste the code and then start using it in formulas.
Dim f_min As Double, f_max As Double, f As Double
Dim g_lower As Double, g_middle As Double
Dim n As Integer
f_min = (2.51 / Re) ^ 2 * (1 - Rr / 3.7) ^ (-2)
If (Rr = 0) Then
f_max = ((Log(10) + 1) / (2 * Log(Re / 2.51))) ^ 2
Else
f_max = ((1 + 2 * 3.7 * 2.51 / (Rr * Re * Log(10))) / (2 * Log(3.7 / Rr) / Log(10))) ^ 2
End If
n = 0
Do
n = n + 1
f = (f_min + f_max) / 2
g_middle = -2 / Log(10) * Log(Rr / 3.7 + 2.51 / Re * f ^ (-1 / 2)) - f ^ (-1 / 2)
g_lower = -2 / Log(10) * Log(Rr / 3.7 + 2.51 / Re * f_min ^ (-1 / 2)) - f_min ^ (-1 / 2)
If (g_middle = 0 Or (f_max - f_min) / 2 < 0.000001) Then
Exit Do
Else
If (g_middle * g_lower > 0) Then
f_min = f
Else
f_max = f
End If
End If
Colebrook = f
If (n > 100) Then
Colebrook = CVErr(xlErrNA)
Exit Do
End If
Loop
End Function
The convergence criteria is the 0.000001 number within the code itself, so you can just change that yourself if desired. I notice that values tend to only be accurate to 10^(-16). In fact, for large Reynolds numbers, the expression for g(fmax) sometimes gives negative values, even though it should only give positive values (this is the entire point of having fmax). But it turns out the cause is that fmax functionally is the answer, and any deviation from zero is due to floating point errors. The root finding routine still gives an answer with Rr up to 3, even though that is mostly unphysical.
For some results, I plugged in a value of Rr of 0.015 and then looked at the output for different Reynolds numbers. I obtained the following, with fmin and fmax from the equations I gave and then the f value from the function. You can see how the upper and lower bounds I reported always bound the solution.
In conclusion, I've shown here a way to get bounds for root finding on the Colebrook equation from the mathematics itself. It also follows that you should be able to count on this method for the full range of inputs, even unphysical ones. For the case of unphysical roughness in laminar regions it does fail, but that's good, because the equation doesn't have a solution there. In all other cases, even when the values are really really close to the unphysical region, this method still works. That's good, because we like things that we can count on working.
Even though it seems like a trivial problem, it's not so easy to get something that won't unexpectedly break, and that's my intention with this approach.
Quicklatex
I just noticed the disclaimer on the website: http://quicklatex.com/
*QuickLaTeX.com is a linkware, free for personal and non-commercial use in exchange to backlink.I have been using this site extensively. It's where basically all of my equations are rendered.
No comments:
Post a Comment