% pipeflowOC.m
% Solution to pipeflow problem using orthogonal collocation method.
% Options include planar, cylindrical or annulus,
% two viscosities (for two immiscible fluids),
% viscosity (Pa s) a function of position and shear rate (for non-Newtonian fluids),
% delta-p/L (Pa/m), total width (H, planar); radius (R, pipe);
% radius (R) and inner radius (R-0, annulus);
% location of interface between two fluids (Hi, planar; Ri, pipe and annulus);
% 0 < Hi < H for planar; 0 < Ri < R for pipe; R-0 < Ri < R for annulus;
% number of grid points in first region (numbered 1 to N),
% second region (numbered N to N+M); total of N+M grid points, interface at N.
% You can solve for one fluid only by using mu1 = mu2 and Hi or Ri 0;
% the orthogonal collocation matrices.
% Planar.m and planar_interp.m must be available.
% Output is velocity and viscosity versus grid points, average velocity
% and flow rate in each region (m^3/s); for planar geometry the flow rate
% is in the flow rate per unit width (into the page), (m^2/s).

global xoc woc aoc boc qinv
% Data input
linear = 1 % linear = 1 if linear, 0 if not
geometry = 1 % planar
%geometry = 2 % pipe
%geometry = 3 % annulus
mu1 = 0.001 % Pa s
mu2 = 0.01 % Pa s
dpL = 1.e3 % Pa/m
if geometry == 1
radius = 0.01 % m
innerradius = 0.
elseif geometry == 2
radius = 0.01 % m
innerradius = 0.
else
radius = 0.01 % m
innerradius = 0.5*radius % m
end
interfaceradius = 0.3*(radius+innerradius) % m
% number of total collocation points in first region
ngrid = 3
% number of total collocation points in second region
mgrid = 3
ntotal = ngrid + mgrid - 1
dr1 = interfaceradius - innerradius
dr2 = radius - interfaceradius
% Call orthogonal collocation matrices
planar(ngrid,1)
xocn=xoc; wocn=woc; aocn=aoc; bocn=boc; qinvn=qinv;
planar(mgrid,1)
xocm=xoc; wocm=woc; aocm=aoc; bocm=boc; qinvm=qinv;

% Need to iterate for nonlinear problem
tol = 1e-12;
% set up initial guess
for i=1:ngrid
x(i) = innerradius + (interfaceradius-innerradius)*xocn(i);
u1(i) = 1;
end
for i=ngrid+1:ntotal
x(i) = interfaceradius + (radius-interfaceradius)*xocm(i-ngrid+1);
u1(i) = 1;
end

%begin iteration
change = 1.;
iter=0;
while change>tol
iter=iter+1;
sum = 0.;
% Set up the matrices
for i=1:ngrid
al(i) = mu1;
end
for i=ngrid+1:ntotal
al(i) = mu2;
end
% set the matrices to zero
for i=1:ntotal
f(i) = 0.;
for j=1:ntotal
aa(i,j) = 0.;
end
end
% set the matrices for the differential equation in the first region
for i=2:ngrid-1
for j=1:ngrid
sum = 0;
for n=1:ngrid
sum = sum + aocn(i,n)*al(n)*aocn(n,j);
end
aa(i,j) = sum;
end
% if geometry > 1
% aa(i,j) = sum + (dr1*al(i)/x(i))*aocn(i,j);
% aa(i,i-1) = aa(i,i-1) - dr1*al(i)/x(i);
% aa(i,i+1) = aa(i,i+1) + dr1*al(i)/x(i);
% end
f(i) = - dr1*dr1*dpL;
end
% need to reset al(ngrid) to the right-hand side value
al(ngrid) = mu2
% set the matrices for the differential equation in the second region
for i=ngrid+1:ntotal-1
for j=ngrid:ntotal
sum = 0;
for n=1:mgrid
sum = sum + aocm(i-ngrid+1,n)*al(n+ngrid-1)*aocm(n,j-ngrid+1);
end
aa(i,j) = sum;
end
% if geometry > 1
% aa(i,j) = sum + (dr2*al(i)/x(i))*aocm(i-ngrid+1,j-ngrid+1);
% aa(i,i-1) = aa(i,i-1) - dr2*al(i)/x(i);
% aa(i,i+1) = aa(i,i+1) + dr2*al(i)/x(i);
% end
f(i) = - dr2*dr2*dpL;
end
% The first and last row of the matrices are special.
f(1) = 0;
f(ntotal) = 0;
% for j=1:ngrid
% aa(1,j) = aocn(1,j);
% end
aa(1,1) = 1.
aa(ntotal,ntotal) = 1;
% The compatibility condition
q = (mu2/mu1)*dr1/dr2
for j=1:ngrid-1
aa(ngrid,j) = aocn(ngrid,j);
end
aa(ngrid,ngrid) = aocn(ngrid,ngrid) - q*aocm(1,1);
for j=ngrid+1:ntotal
aa(ngrid,j) = - q*aocm(1,j-ngrid+1);
end
f(ngrid) = 0;

% Solve one iteration.
u2 = aa\f';
% Calculate the criterion to stop the iterations.
sum=0.;
for i=1:ntotal
sum = sum + abs(u1(i) - u2(i));
u1(i) = u2(i);
end
if linear==1,break,end
change = sum
end

% Plot the solution.
% Split the solution into two regions and find the coefficients
for i=1:ngrid
vel1(i) = u1(i);
x1(i) = x(i);
end
dd1 = qinvn*vel1';
xp1 = 0:dr1/20:dr1; %plotting points
uinterp1 = planar_interp(ngrid,dd1,xp1,0.,dr1);
for i=ngrid:ntotal
vel2(i-ngrid+1) = u1(i);
x2(i-ngrid+1) = x(i);
end
dd2=qinvm*vel2';
xp2 = 0:dr2/20:dr2;
xp2 = dr1 + xp2; %plotting points
uinterp2 = planar_interp(mgrid,dd2,xp2,dr1,dr1+dr2);
plot(x,u1,'*b',xp1,uinterp1,'-b',xp2,uinterp2,'-b')
xlabel('x')
ylabel('u')
disp(u1)

% calculate the area
sum = 0.;
sum1 = 0.;
for i=1:ngrid
sum = sum + wocn(i);
sum1 = sum1 + wocn(i)*u1(i);
end
areaone = (interfaceradius - innerradius)*sum
flow1 = (interfaceradius - innerradius)*sum1
sum = 0.;
sum1 = 0.;
for i=ngrid:ntotal
sum = sum + wocm(i-ngrid+1);
sum1 = sum1 + wocm(i-ngrid+1)*u1(i);
end
areatwo = (radius - interfaceradius)*sum
flow2 = (radius - interfaceradius)*sum1
avgvel1 = flow1/areaone
avgvel2 = flow2/areatwo
fractionflow1 = flow1/(flow1+flow2)