Skip to content

Instantly share code, notes, and snippets.

@alexlib

alexlib/colebrook.ipynb

Created November 18, 2012 21:51
Show Gist options
  • Save alexlib/4107731 to your computer and use it in GitHub Desktop.
Save alexlib/4107731 to your computer and use it in GitHub Desktop.
Download ZIP
Translation of the Matlab function colebrook.m to IPython notebook
Display the source blob
Display the rendered blob
Raw
colebrook.ipynb
{
"metadata": {
"name": "colebrook"
},
"nbformat": 2,
"worksheets": [
{
"cells": [
{
"cell_type": "markdown",
"source": [
"# Translation of Matlab function",
"### Fast and accurate computation of the Darcy-Weisbach friction factor F according to the Colebrook-White equation:",
"$$\\frac{1}{\\sqrt{F}} = -2 \\log_{10} \\left\\[ \\frac{K}{3.7} + \\frac{2.51}{R \\sqrt{F}} \\right\\]$$",
"",
"<html>",
"<code>",
"function F = colebrook(R,K)",
"% F = COLEBROOK(R,K) fast, accurate and robust computation of the ",
"% Darcy-Weisbach friction factor F according to the Colebrook equation:",
"% - -",
"% 1 | K 2.51 |",
"% --------- = -2 * Log_10 | ----- + ------------- |",
"% sqrt(F) | 3.7 R * sqrt(F) |",
"% - -",
"% INPUT:",
"% R : Reynolds' number (should be >= 2300).",
"% K : Equivalent sand roughness height divided by the hydraulic ",
"% diameter (default K=0).",
"%",
"% OUTPUT:",
"% F : Friction factor.",
"%",
"% FORMAT:",
"% R, K and F are either scalars or compatible arrays.",
"%",
"% ACCURACY:",
"% Around machine precision forall R > 3 and forall 0 <= K, ",
"% i.e. forall values of physical interest. ",
"%",
"% EXAMPLE: F = colebrook([3e3,7e5,1e100],0.01)",
"%",
"% Edit the m-file for more details.",
"",
"% Method: Quartic iterations.",
"% Reference: http://arxiv.org/ ",
"% Read this reference to understand the method and to modify the code.",
"",
"% Author: D. Clamond, 2008-09-16. ",
"",
"% Check for errors.",
"if any(R(:)<2300) == 1, ",
" warning('The Colebrook equation is valid for Reynolds'' numbers >= 2300.'); ",
"end,",
"if nargin == 1 | isempty(K) == 1, ",
" K = 0;",
"end,",
"if any(K(:)<0) == 1, ",
" warning('The relative sand roughness must be non-negative.'); ",
"end,",
"",
"% Initialization.",
"X1 = K .* R * 0.123968186335417556; % X1 <- K * R * log(10) / 18.574.",
"X2 = log(R) - 0.779397488455682028; % X2 <- log( R * log(10) / 5.02 ); ",
"X1+X2; ",
"% Initial guess. ",
"F = X2 - 0.2; ",
"",
"% First iteration.",
"E = ( log(X1+F) - 0.2 ) ./ ( 1 + X1 + F );",
"F = F - (1+X1+F+0.5*E) .* E .* (X1+F) ./ (1+X1+F+E .* (1+E/3));",
"",
"% Second iteration (remove the next two lines for moderate accuracy).",
"E = ( log(X1+F) + F - X2 ) ./ ( 1 + X1 + F );",
"F = F - (1+X1+F+0.5*E) .* E .* (X1+F) ./ (1+X1+F+E .* (1+E/3));",
"",
"% Finalized solution.",
"F = 1.151292546497022842 ./ F; % F <- 0.5 * log(10) / F;",
"F = F .* F; % F <- Friction factor.",
"</code>",
"</html>"
]
},
{
"cell_type": "markdown",
"source": []
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"from math import log",
"",
"def colebrook(R,K):",
" '''",
"% function F = colebrook(R,K)",
"% F = COLEBROOK(R,K) fast, accurate and robust computation of the ",
"% Darcy-Weisbach friction factor F according to the Colebrook equation:",
"% - -",
"% 1 | K 2.51 |",
"% --------- = -2 * Log_10 | ----- + ------------- |",
"% sqrt(F) | 3.7 R * sqrt(F) |",
"% - -",
"% INPUT:",
"% R : Reynolds' number (should be >= 2300).",
"% K : Equivalent sand roughness height divided by the hydraulic ",
"% diameter (default K=0).",
"%",
"% OUTPUT:",
"% F : Friction factor.",
"%",
"% FORMAT:",
"% R, K and F are either scalars or compatible arrays.",
"%",
"% ACCURACY:",
"% Around machine precision forall R > 3 and forall 0 <= K, ",
"% i.e. forall values of physical interest. ",
"%",
"% EXAMPLE: F = colebrook([3e3,7e5,1e100],0.01)",
"%",
"% Edit the m-file for more details.",
"",
"",
"% Method: Quartic iterations.",
"% Reference: http://arxiv.org/ ",
"% Read this reference to understand the method and to modify the code.",
"",
"",
"% Author: D. Clamond, 2008-09-16. ",
"",
"",
"% Check for errors.",
"if any(R(:)<2300) == 1, ",
" warning('The Colebrook equation is valid for Reynolds'' numbers >= 2300.'); ",
"",
"end,",
"if nargin == 1 | isempty(K) == 1, ",
"",
" K = 0;",
"end,",
"if any(K(:)<0) == 1, ",
" warning('The relative sand roughness must be non-negative.'); ",
"end,",
"",
"",
"% Initialization.",
"X1 = K .* R * 0.123968186335417556; % X1 <- K * R * log(10) / 18.574.",
"X2 = log(R) - 0.779397488455682028; % X2 <- log( R * log(10) / 5.02 ); ",
"",
"X1+X2; ",
"% Initial guess. ",
"",
"F = X2 - 0.2; ",
"",
"",
"% First iteration.",
"E = ( log(X1+F) - 0.2 ) ./ ( 1 + X1 + F );",
"F = F - (1+X1+F+0.5*E) .* E .(X1+F) ./ (1+X1+F+E.(1+E/3));",
"",
"",
"% Second iteration (remove the next two lines for moderate accuracy).",
"E = ( log(X1+F) + F - X2 ) ./ ( 1 + X1 + F );",
"F = F - (1+X1+F+0.5*E) .* E .(X1+F) ./ (1+X1+F+E.(1+E/3));",
"",
"",
"% Finalized solution.",
"F = 1.151292546497022842 ./ F; % F <- 0.5 * log(10) / F;",
"F = F .* F; % F <- Friction factor.",
"'''",
" ",
" if R < 2300: ",
" raise ValueError('The Colebrook equation is valid for Reynolds'' numbers >= 2300.') ",
"",
" if K == []:",
" K = 0",
" elif K < 0:",
" raise ValueError('The relative sand roughness must be non-negative.') ",
"",
"",
" X1 = K * R * 0.123968186335417556 # X1 <- K * R * log(10) / 18.574.",
" X2 = log(R) - 0.779397488455682028 # X2 <- log( R * log(10) / 5.02 ); ",
"",
"",
" F = X2 - 0.2; ",
"",
"",
" # First iteration.",
" E = ( log(X1+F) - 0.2 ) / ( 1 + X1 + F )",
" F = F - (1+X1+F+0.5*E) * E * (X1+F)/(1+X1+F+E*(1+E/3))",
"",
" # Second iteration (remove the next two lines for moderate accuracy).",
" E = ( log(X1+F) + F - X2 ) / ( 1 + X1 + F )",
" F = F - (1+X1+F+0.5*E) * E * (X1+F) / (1+X1+F+E*(1+E/3))",
"",
"",
" # Finalized solution.",
" F = 1.151292546497022842 / F # F <- 0.5 * log(10) / F;",
" F = F * F # F <- Friction factor.",
" ",
" return F"
],
"language": "python",
"outputs": [],
"prompt_number": 26
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"for r in [3e3,7e5,1e100]:",
" print colebrook(r,0.01)"
],
"language": "python",
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"0.0518683608506",
"0.0379908240717",
"0.0379037118924"
]
}
],
"prompt_number": 28
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"help(colebrook)"
],
"language": "python",
"outputs": [
{
"output_type": "stream",
"stream": "stdout",
"text": [
"Help on function colebrook in module __main__:",
"",
"colebrook(R, K)",
" % function F = colebrook(R,K)",
" % F = COLEBROOK(R,K) fast, accurate and robust computation of the ",
" % Darcy-Weisbach friction factor F according to the Colebrook equation:",
" get_ipython().magic(u'- -')",
" % 1 | K 2.51 |",
" % --------- = -2 * Log_10 | ----- + ------------- |",
" % sqrt(F) | 3.7 R * sqrt(F) |",
" % - -",
" % INPUT:",
" get_ipython().magic(u\"R : Reynolds' number (should be >= 2300).\")",
" % K : Equivalent sand roughness height divided by the hydraulic ",
" % diameter (default K=0).",
" %",
" % OUTPUT:",
" get_ipython().magic(u'F : Friction factor.')",
" %",
" % FORMAT:",
" get_ipython().magic(u'R , K and F are either scalars or compatible arrays.')",
" %",
" % ACCURACY:",
" get_ipython().magic(u'Around machine precision forall R > 3 and forall 0 <= K,')",
" % i.e. forall values of physical interest. ",
" %",
" % EXAMPLE: F = colebrook([3e3,7e5,1e100],0.01)",
" %",
" % Edit the m-file for more details.",
" ",
" ",
" % Method: Quartic iterations.",
" % Reference: http://arxiv.org/ ",
" % Read this reference to understand the method and to modify the code.",
" ",
" ",
" % Author: D. Clamond, 2008-09-16. ",
" ",
" ",
" % Check for errors.",
" if any(R(:)<2300) == 1, ",
" warning('The Colebrook equation is valid for Reynolds'' numbers >= 2300.'); ",
" ",
" end,",
" if nargin == 1 | isempty(K) == 1, ",
" ",
" K = 0;",
" end,",
" if any(K(:)<0) == 1, ",
" warning('The relative sand roughness must be non-negative.'); ",
" end,",
" ",
" ",
" % Initialization.",
" X1 = K .* R * 0.123968186335417556; % X1 <- K * R * log(10) / 18.574.",
" X2 = log(R) - 0.779397488455682028; % X2 <- log( R * log(10) / 5.02 ); ",
" ",
" X1+X2; ",
" % Initial guess. ",
" ",
" F = X2 - 0.2; ",
" ",
" ",
" % First iteration.",
" E = ( log(X1+F) - 0.2 ) ./ ( 1 + X1 + F );",
" F = F - (1+X1+F+0.5*E) .* E .(X1+F) ./ (1+X1+F+E.(1+E/3));",
" ",
" ",
" % Second iteration (remove the next two lines for moderate accuracy).",
" E = ( log(X1+F) + F - X2 ) ./ ( 1 + X1 + F );",
" F = F - (1+X1+F+0.5*E) .* E .(X1+F) ./ (1+X1+F+E.(1+E/3));",
" ",
" ",
" % Finalized solution.",
" F = 1.151292546497022842 ./ F; % F <- 0.5 * log(10) / F;",
" F = F .* F; % F <- Friction factor.",
""
]
}
],
"prompt_number": 30
},
{
"cell_type": "code",
"collapsed": true,
"input": [],
"language": "python",
"outputs": []
}
]
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment