% DOL 10/13/09
% Octave (Matlab) script to solve Type III pipe flow proplem.
% Find Diameter, given head loss.  
% Use a system demand curve.

clc;
clear;

%---------- Given these:

mu = 1.985E-5;      % kg/m*s
rho = 1.145;        % kg/m3
L   = 150;          % m
eps = 0.0;          % m
Q   = 0.35;         % m3/s
hLwant = 20;        % m, desired head loss

%---------

g = 9.80665;        % m/s2
D = linspace(0.01,1, 1000);
D = D';

v = Q./(pi/4*D.^2);
Re = rho/mu*D.*v;

%--------- Compute f

% Initial guess from Haaland equation; x = sqrt(f)

x = 1.0./(-1.8*log10(6.9./Re + (eps./D/3.7).^1.11));  % x is sqrt(f)

% Solve Colbrook equation with Newton's method, 5 iterations (plenty)
% Newtons method is x = x0 - F0/(dF/dx)_0
% Here, F(x) = 0; with F(x) = 1/x+2*log10(eps/3.7/D + 2.51/Re/x)
% This is just the Colbrook equation rearranged to equal zero

for i=1:5       
    x = x + (1./x + 2*log10(eps./D/3.7 + 2.51./Re./x)) ./ ...
            (1./x./x + 2*2.51./Re./x./x./log10(eps./D/3.7+2.51./Re./x));
end

% Recover f

f = x.^2;   

%---------- Compute head loss

hL = f*L./D .*v.^2/2/g;

%---------- Plot the results

semilogy(D,hL);
xlabel('D (m)');
ylabel('hL (m)');

%---------- interpolate D for desired head loss

Dwant = interp1(hL,D,hLwant)

