% pipeflowNR.m
% Solution to pipeflow problem using finite difference method, version3
% 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 finite difference method uses a continuity equation at grid point N
% rather than the differential equation.
% 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).

% 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.010 % 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.5*(radius+innerradius)% m
ngrid = 13
mgrid = 12
ntotal = ngrid+mgrid
dr1 = (interfaceradius - innerradius)/(ngrid-1)
dr2 = (radius - interfaceradius)/mgrid

% Need to iterate for nonlinear problem
tol = 1e-12;
% set up initial guess
for i=1:ngrid
u1(i) = innerradius + dr1*(i-1);
x(i) = innerradius + dr1*(i-1);
end
for i=ngrid+1:ngrid+mgrid
u1(i) = interfaceradius + dr2*i;
x(i) = interfaceradius + dr2*(i-ngrid);
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 i=2:ngrid-1
aa(i,i) = -(al(i+1)+2*al(i)+al(i-1));
aa(i,i-1) = al(i)+al(i-1);
aa(i,i+1) = al(i)+al(i+1);
if geometry > 1
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) = - 2*dr1*dr1*dpL;
end
% need to reset al(ngrid) to the right-hand side value
al(ngrid) = mu2
for i=ngrid+1:ntotal-1
aa(i,i) = -(al(i+1)+2*al(i)+al(i-1));
aa(i,i-1) = al(i)+al(i-1);
aa(i,i+1) = al(i)+al(i+1);
if geometry > 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) = - 2*dr2*dr2*dpL;
end
% The first and last row of the matrices are special.
f(1) = 0;
f(ntotal) = 0;
aa(1,1) = 1;
aa(1,2) = 0;
aa(ntotal,ntotal-1) = 0;
aa(ntotal,ntotal) = 1;
% The compatibility condition
p = mu1
q = mu2*dr1/dr2
aa(ngrid,ngrid-1) = -4*mu1;
aa(ngrid,ngrid+1) = -4*q;
pfactor = -p/aa(ngrid-1,ngrid-2)
qfactor = -q/aa(ngrid+1,ngrid+2)
aa(ngrid,ngrid-1) = aa(ngrid,ngrid-1) + pfactor*aa(ngrid-1,ngrid-1);
aa(ngrid,ngrid) = 3*(p+q)+pfactor*aa(ngrid-1,ngrid-2)+qfactor*aa(ngrid+1,ngrid);
aa(ngrid,ngrid+1) = aa(ngrid,ngrid+1) + qfactor*aa(ngrid+1,ngrid+1);
f(ngrid) = pfactor*f(ngrid-1) + qfactor*f(ngrid+1);
% 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 % do only one iteration
change = sum
end

% Plot the solution.
plot(x,u1,'-')
xlabel('x')
ylabel('u')
u1
% put the two regions into different vectors for integration
for i=1:ngrid
u11(i) = u1(i);
x1(i) = x(i);
area1(i) = 1.;
end
for i=ngrid:ntotal
u12(i-ngrid+1) = u1(i);
x2(i-ngrid+1) = x(i);
area2(i-ngrid+1) = 1.;
end
% calculate the area
areaone = trapz(x1,area1)
areatwo = trapz(x2,area2)
flow1 = trapz(x1,u11)
flow2 = trapz(x2,u12)
avgvel1 = flow1/areaone
avgvel2 = flow2/areatwo
fractionflow1 = flow1/(flow1+flow2)