function NewtonFluidFlow
% 水平管中不可壓縮牛頓流體的層流流動
% 求解動量守衡方程以獲得剪切力和速度分布
% Laminar Flow in a Horizontal Pipe (Newtonian Fluid)
%
%   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

global deltaP L mu R
deltaP = 500;
L = 10;
mu = 8.937e-4;
R = 0.009295;
a = 0;
b = R;

% 求解ODE-BVP問題
% initialize of solution with a guess of y1(r)=0,y2(r)=0,y3(r)=0  
solinit = bvpinit(linspace(a,b,100),[0 0 0]);
sol = bvp4c(@ODEs,@BCfun,solinit);

% 分析解
r = sol.x;
TauAnal = (deltaP/(2*L))*r;
uAnal = (deltaP*R^2/(4*mu*L))*(1-(r/R).^2);
umAnal = deltaP*R^2/(8*mu*L);

% 繪制圖
% 剪應力分布
tau = sol.y(1,2:end)./sol.x(2:end);
tau = [0 tau];
plot(r,tau,'b-',r,TauAnal,'r-.')
xlabel('管徑位置r，m')
ylabel('剪應力，kg/(m s^2)')
legend('數值解','分析解')
figure
% 剪切速率分布
plot(r,sol.y(2,:),'b-',r,uAnal,'r-.')
xlabel('管徑位置r，m')
ylabel('剪切速率，m/s')
legend('數值解','分析解')

% 計算催化劑的有效因子
fprintf('\t平均剪切速率為: um = %.4f',sol.y(3,end))

% ------------------------------------------------------------------
function dydr = ODEs(r,y)
global deltaP L mu R
rTau = y(1);
u = y(2);
um = y(3);
drTaudr = deltaP*r/L; 
if r > 0             % 避免除數為零
    Tau = rTau/r;
else
    Tau = 0;
end
dudr = -Tau/mu;
dumdr = u*2*r/R^2;

dydr = [drTaudr; dudr; dumdr];

% ------------------------------------------------------------------
function bc = BCfun(ya,yb)
bc = [ya(1); yb(2); ya(3)];