Calculate steady-state boundary shear stress in pressurised water pipes

shear_stress.m

function tau_a = shear_stress(D, Q, k_s, T, den)
    %{ Hydraulic shear stress at pipe wall (in Pa).
    %
    % Params:
    % D -- internal diameter (m)
    % Q -- flow (m^3s^-1)
    % k_s -- roughness height (m)
    % T -- temperature; defaults to 10degC)
    % den -- density defaults to 1000kg/m^3
    % 
    % Will Furnass (wrfurnass1@sheffield.ac.uk)
    % Pennine Water Group, University of Sheffield
    % November 2014
    %
    % Ported from Python (see https://github.com/willfurnass/pyhyd)
    % License: GPL v3 (www.gnu.org/copyleft/gpl.html)
    %}
    if not(exist('T', 'var'))
        T = 10.0;
    end
    if not(exist('den', 'var'))
        den = 1000.0;
    end
    if any(den <= 0)
        throw(MException('shear_stress:invalid_argument', ...
                         'Non-positive density'));
    end
    if any(D <= 0)
        throw(MException('shear_stress:invalid_argument', ...
                         'Non-positive internal pipe diam.'));
    end
    g = 9.81;
    tau_a = den * g * (D / 4.0) * ...
            hyd_grad(D, Q, k_s, T, den);
end

function A_xc = x_sec_area(D)
    %{ Cross-sectional area of a uniform pipe
    % 
    % Params:
    % D -- internal diameter 
    %}
    if any(D <= 0)
        throw(MException('shear_stress:invalid_argument', ... 
                         'Non-positive internal pipe diam.'));
    end
    A_xc =  pi * power(D, 2) / 4.0;
end


function mu = dyn_visc(T)
    %{ Dynamic viscosity of water (N.s.m^-2)
    %
    % Params:
    % T -- temperature (degC) (defaults to 10degC)
    %
    % See http://en.wikipedia.org/wiki/Viscosity#Viscosity_of_water for eq.
    %}
    if not(exist('T', 'var'))
        T = 10.0;
    end

    if any(T < 0 | T > 100)
        err_msg = strcat('Cannot calc dynamic viscosity: ', ...
                         'temperature outside range [0,100].');
        throw(MException('shear_stress:invalid_argument', err_msg)); 
    end
    A = 2.414e-5; % Pa.s
    B = 247.8;  % K
    C = 140.0;  % K
    mu = A * power(10, (B / (T + 273.15 - C)));
end


function Re = reynolds(D, Q, T, den)
    %{Reynolds number
    %
    % Params:
    % D -- internal diameter (m)
    % Q -- flow (m^3s^-1)
    % T -- temperature; defaults to 10degC
    % den -- density; defaults to 1000kg/m^3
    %}
    if not(exist('T', 'var'))
        T = 10.0;
    end
    if not(exist('den', 'var'))
        den = 1000.0;
    end
    if any(den <= 0)
        throw(MException('shear_stress:invalid_argument', ...
            'Non-positive fluid density.'));
    end
    Re = ((abs(Q) / x_sec_area(D)) * D) / (dyn_visc(T) / den);
end


function f = friction_factor(D, Q, k_s, T, den)
    %{ Darcy-Weisbach friction factor.
    % 
    % Params:
    % D -- internal diameter (m)
    % Q -- flow (m^3s^-1)
    % k_s -- roughness height (m)
    % T -- temperature; defaults to 10degC)
    % den -- density defaults to 1000kg/m^3
    %
    % Laminar flow:      Hagen-Poiseuille formula
    % Transitional flow: cubic interpolation from Moody Diagram
    %                    for transition region as per the EPANET2 manual
    %                    (in turn taken from Dunlop (1991))
    % Turbulent flow: Swamee-Jain approximation of implicit
    %                 Colebrook-White equation
    %}
    if not(exist('T', 'var'))
        T = 10.0;
    end
    if not(exist('den', 'var'))
        den = 1000.0;
    end
    if any(k_s < 0)
        throw(MException('shear_stress:invalid_argument', ...
            'Negative pipe roughness.'));
    end
    Re = reynolds(D, Q, T, den);
    if Re == 0 
        f = 0 * Re;
    elseif Re < 2000
        f = 64 / Re;
    elseif (2000 <= Re) && (Re < 4000)
        y3 = -0.86859 * log((k_s / (3.7 * D)) + (5.74 / (4000^0.9)));
        y2 = (k_s / (3.7 * D)) + (5.74 / power(Re, 0.9));
        fa = power(y3, -2);
        fb = fa * (2 - (0.00514215 / (y2*y3)));
        r = Re / 2000.;
        x4 = r * (0.032 - (3. * fa) + (0.5 * fb));
        x3 = -0.128 + (13. * fa) - (2.0 * fb);
        x2 = 0.128 - (17. * fa) + (2.5 * fb);
        x1 = (7 * fa) - fb;
        f = x1 + r * (x2 + r * (x3 + x4));
    elseif Re >= 4000
        %if warn
        %    if k_s < 4e-4 || k_s > 5e-2:
        %    err_msg = sprintf(strcat('Swamee Jain approx to Colebrook White not ', ...
        %                             'valid for turb flow as k_s=%fm ', ...
        %                             '(outside range [0.004,0.05])'), k_s);
        %        throw(MException('shear_stress:invalid_argument',
        %        err_msg));
        %    end
        %    if Re > 1e7
        %        err_msg = sprintf(strcat('Swamee Jain approx to Colebrook White not ', ...
        %                                 'valid for turb flow as Re=%f ', ...
        %                                 'greater than 10,000,000'), Re);
        %        throw(MException('shear_stress:invalid_argument', err_msg));
        %    end
        %end
        f = 0.25 / power((log10((k_s / (3.7 * D)) + ...
                (5.74 / power(Re, 0.9)))), 2);
    end
end


function S_0 = hyd_grad(D, Q, k_s, T, den)
    %{ Headloss per unit length of pipe (in m), using approx. to Colebrook White eq.
    %
    % Params:
    % D -- internal diameter (m)
    % Q -- flow (m^3s^-1)
    % k_s -- roughness height (m)
    % T -- temperature; defaults to 10degC)
    % den -- density defaults to 1000kg/m^3
    %}
    if not(exist('T', 'var'))
        T = 10.0;
    end
    if not(exist('den', 'var'))
        den = 1000.0;
    end
    if any(den <= 0)
        throw(MException('shear_stress:invalid_argument', ...
            'Non-positive fluid density.'));
    end
    if any(D <= 0)
        throw(MException('shear_stress:invalid_argument', ...
            'Non-positive internal pipe diam.'));
    end
    f = friction_factor(D, Q, k_s, T, den);
    vel_sq = power((Q / x_sec_area(D)), 2);
    g = 9.81;
    S_0 = (f * vel_sq) / (D * 2 * g);
end