function PipeDiameter
% 估算非牛頓型流體流動所需的管道直徑
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/05/30 $

clear all
clc

% 一水平管道內冪定律流體的流動
rho = 961;
G = 6.67;
L = 10;
eps = 5e-6;
deltaP = 15e3;
k = 1.8;
n = 0.64;
K = 1.48;
n1 = n;
K1 = K*((3*n+1)/(4*n))^n;

D1 = (32*K1*L*8^(n1-1)*(4*G/pi/rho)^n1)^(1/(3*n1+1));
D2 = 2*D1;
delta = 1e-4;
while abs((D2 - D1)/D2) > delta
    u = (G/rho)/(pi*D1^2/4);
    Regen = D1^n1*u^(2-n1)*rho/(8^(n1-1)*K1);
    if Regen <= 2100
        if k == 0
            break
        else
            f = 16/Regen;
        end
    end
    if Regen > 2100
        f = FrictFactor(Regen,n1);
    end
    Le = k*D1/(4*f);
    D1 = D2;
    D2 = (2*f*(L+Le)/(rho*deltaP)*(4*G/pi)^2)^0.2;
end
fprintf('\t管道直徑為D = %.4f %s\n',D2,'m')
fprintf('\t摩擦因子為f = %.4f %s',f,'m')

% ------------------------------------------------------------------
function f = FrictFactor(Regen,n1)
f0 = 16/Regen;      %4.53
f = fzero(@fFunc,f0,[],Regen,n1)

% ------------------------------------------------------------------
function F = fFunc(f)
F = 1/sqrt(f) - 4.0*log(Regen*f^(1-n1/2))/n1^0.75 + 0.4/n1^1.2;%4.54