% Time-dependent Hagen-Poiseuille flow simulation
% EPS 109 Fall 2011, Final Project
% Prof. Burkhard Militzer, GSI Laura Nielson
% Aaron Tran

% Start from the Navier-Stokes equations
% Incompressible, Newtonian, axisymmetric flow
% In a cylinder with linear pressure gradient
% Velocity is only directed in the forward (z) direction
% Thus U_r = U_phi = 0 everywhere - there is no radial or tangential flow

% The continuity equation becomes $\partial U_z / \partial z = 0$
% Thereby eliminating the need to solve over a 2-D grid
% We can solve a 1-D PDE:
% \rho \frac{\partial u_z}{\partial t} = -\alpha + \mu (\frac{1}{r} \frac{\partial u_z}{\partial r} + \frac{\partial^2 u_z}{\partial r^2})

clear
clf

% Constants
mu = 1; % Dynamic viscosity (Pa s)
rho = 1000; % Fluid density (kg/m^3)
alpha = -1; % Pressure gradient (Pa/m) satisfying P(z) = P(0) + alpha*z
R = 1; % Pipe radius (m)

% Simulation parameters
dt = 0.05; % Time step (s)
dr = 0.05; % Spatial step over the radius (m)
N = R/dr; % Number of intervals in grid (please set dr, R so that this is an integer)
iter = 30000; % Number of iterations

frameint = 150; % Number of iterations between frame captures (for movie)

% Initial velocity profile, v_z(r)
% We only do computations for half of the pipe and double it in the visualization step
% We move from the center of the pipe to the outside, so v(1) is at the exact center
v = zeros(N+1,1);
v(1:N) = 0; % Initial velocity distribution
for i = 1:N % Funky initial velocity distribution
    v(i) = 0.2*cos(pi*i/N);
end
v(N+1) = 0.1; % Edge velocity, with no-slip (Dirichlet) boundary
vnew = zeros(N+1,1);
vnew(N+1) = v(N+1);

% Both 1st and 2nd derivative stencils are centered
% This is prone to instability and could be much improved
f = 0; % Image numbering
for i = 1:iter
    
	% Visualization step - done first to get initial conditions
	if mod(i-1,frameint) == 0
        plot([flipud(v);v(2:end)],(0:dr:2*R));
    	hold on
        line([-10 10],[2*R 2*R],'LineWidth',3,'Color','k') % [x1 x2] [y1 y2]
        line([-10 10],[0 0],'LineWidth',3,'Color','k')
        hold off
        
        title(sprintf('Time evolving Poiseuille flow, t = %0.1f seconds',(i-1)*dt));
        xlabel('Velocity (m/s)');
        ylabel('Cross sectional position of pipe (m)');
        axis([-0.1 0.5 -R/2 5/2*R]); % Hardcoded for v axis, adjust as needed
        
        fprintf('Iteration %d of %d\n',i,iter);
        drawnow
        
% 		saveas(gcf, ['final_pres_var',sprintf('%03d',f),'.png'])
        f = f + 1;
    end
    
    % Step forward in time
	for j = 1:N % Don't change our boundary
		r_pos = (j-1)*dr;
        if r_pos ~= 0
			vnew(j) = v(j) + dt/rho * (-alpha + mu*((v(j+1)-v(j-1))/(2*dr*r_pos) + (v(j+1)-2*v(j)+v(j-1))/(dr^2) ));
        else % j = 1, we are at center
			vnew(j) = v(j) + dt/rho * (-alpha + mu*(2*v(j+1)-2*v(j))/(dr^2)); % Reflect across center, discard the 1/r term (take the limit)
        end
	end
    v = vnew;
end

% Plot the steady-state solution for comparison
y = (0:dr:2*R);
x = -1/(4*mu) * alpha * (R^2 - (y-R).^2);
hold on
plot(x,y,'r-o')

title(sprintf('Time evolving Poiseuille flow, t = %0.1f seconds',i*dt));
% saveas(gcf, ['final_pres_var',sprintf('%03d',f),'.png'])