%Project 4: Unsteady-state motion-driven flow between two parallel plates
%           with lower plate moving and upper plate held stationary

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

clear all

iF = 11; %iF is the number of nodes in space
L = 1; %nondimensional distance between the plates
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

jF = 100; %jF is the number of nodes in time
sumerr = 1;
while sumerr > 0.01
    jF = jF+100;
    
    %Boundary conditions
    for j = 1:jF
        U(1,j) = 1; %lower plate set in motion
        U(iF,j) = 0; %upper plate held boundary
    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);
        end
    end
    err = U(:,jF)-U(:,jF-100);
    sumerr = err'*err
end    

%Simulating the dimensional velocity profile
B = 0.1; %width of the fluid 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)
V = 2; %velocity of the moving plate in m/sec

for i = 1:iF
    y(i) = (i-1)*is*B;
    for j = 1:jF
        t(j) = (j-1)*js*rho*B*B/mu;
        u(i,j) = V*U(i,j);
    end
    uss(i) = V*(1-y(i)/B); %steady velocity profile
 end


%Simulating the volumetric flow rate
for j = 1:jF
    Q(j) = 0;
    for i = 2:iF
        Q(j) = Q(j) + 0.5*(u(i,j)+u(i-1,j))*is*B;
    end
end
Qss = B*V/2 %steady volumetric flow rate

subplot(2,1,1), plot(y,u,y,uss,'.')
xlabel('Distance from the moving lower plate to the stationary upper plate (in m)')
ylabel('Velocity (in m/s)')
title('Project 4: Motion-driven flow between two parallel plates')
subplot(2,1,2), plot(t/60,Q,t/60,Qss,'.')
xlabel('time (in min)')
ylabel('Volumetric flow (in m^3/s)')
legend('unsteady profile','steady profile','location','best')