%Project 2 CORRECTED: Unsteady-state pressure-driven flow in a circular pipe 

%Explicit method for parabolic PDE is used where c = 1

clear all

iF = 11; %iF is the number of nodes in space
L = 1; %nondimesional radius of the cylinder
is = L/(iF-1); %is the step size in space

% Choosing js (the step size in time)
% r = c*js/(is*is) should be positive and less than 0.5 to ascertain stability
c = 1; %coefficient of the second derivative in the spatial coordinate 
r = 0.4; %r is chosen to be less than 0.5
js = r*(is*is)/c; %js is the step size in time

%Initial condition
for i = 1:iF
    U(i,1) = 0;
end

%Defining the nondimesional length 
for i = 1:iF
    neta(i) = (i-1)*is;
end

jF = 100; %jF is the number of nodes in time
sumerr = 1;
while sumerr > 0.001
    jF = jF+100;
   
    %Boundary condition 1
    for j = 1:jF
        U(1,j) = 0; %no slip at the stationary cylinder wall
    end

    for j = 1:jF-1
        for i = 2:iF-1
            U(i,j+1) = r*U(i-1,j)+(1-2*r)*U(i,j)+r*U(i+1,j)+0.5*(js/is)*(1/(neta(i)-1))*(U(i+1,j)-U(i-1,j))+js*4;
        end
        U(iF,j+1) = U(iF-1,j+1);%first derivative of velocity is zero at the axis(Boundary condition 2)
    end
    err = U(:,jF)-U(:,jF-100);
    sumerr = err'*err
end    

%Simulating the dimensional velocity profile
DelP = 1; %pressure gradiet in kg/(m.s^2)
R = 0.2; %radius of the outer pipe in m
rho = 1000; %density of water at 27 deg C in kg/m^3
mu = 8.90/10000 ; %viscosity of water at 27 deg C in kg/(m.s)

for i=1:iF
    r(i) = (1-neta(i))*R;
    for j = 1:jF
        t(j) = (j-1)*js*rho*R*R/mu;
        u(i,j) = (DelP*R^2/(4*mu))*U(i,j);
    end
    uss(i) = (DelP*R^2/(4*mu))*(1-(r(i)/R)^2);
end

%Simulating the volumetric flow rate
for j = 1:jF
    Q(j) = 0;
    for i = 2:iF
        Q(j) = Q(j) + 2*pi*0.5*(r(i)*u(i,j)+r(i-1)*u(i-1,j))*is*R;
    end
end

Qss = (pi*R^4*DelP/(8*mu)) %steady volumetric flow rate

subplot(2,1,1), plot(r,u,r,uss,'^')
xlabel('Distance from the central axis to the cylinder wall (in m)')
ylabel('Velocity (in m/s)')
title('Project 2: Pressure-driven flow in a circular pipe ')
subplot(2,1,2), plot(t/60,Q,t/60,Qss,'^')
xlabel('Time (in min)')
ylabel('Steady flow rate (in m^3/s)')