% 10.34 - Fall 2006
% Hw Set #8 - Problem #1
% Rob Ashcraft - Oct. 25, 2006


% Heat transfer to a laminar flow in a tube, without neglecting axial
% heat conduction.

%function pset8p1_main

clear; clc;

% Peclet number (Pe) = 2*U*R/a = 2*U*R*rho*Cp/lambda

Re = 100;

T0 = 300;  % K
T1 = 400;  % K
Cp = 4000;  % J/kg-K
lam = 0.6;  % W/m-K = J/m-s-K
mu = 0.5e-3;  % N-s/m^2
rho = 1000;  % kg/m^3

alpha = lam/rho/Cp  % thermal diffusivity (m^2/s)

R = 0.050;  % m

U = Re*mu/(2*R*rho)  % m/s

Pe = 2*U*R/alpha   % Peclet number

tau_cond = R^2/(2*alpha)   % sec

Nr = 21;
Nz = 201;
N_total = Nz*Nr;

Nup = 0.1;
Ndown = 2;
% actual variables
% take X tau's worth of U_centerline
z_initial = -Nup*2*U*tau_cond   %meters, start 2 tau's upstream
z_final = Ndown*2*U*tau_cond   %meters, go 10 tau's downstream

dz = (z_final - z_initial)/(Nz-1);
dr = R/(Nr-1);

% dimensionless variables
zeta_initial = -Nup;   %start 2 tau's upstream
zeta_final = Ndown;   %go 10 tau's downstream

dzeta = (zeta_final - zeta_initial)/(Nz-1);
deta = 1/(Nr-1);

A_sparse = spalloc(N_total,N_total,5*N_total);
b_vec = zeros(Nr*Nz,1);

for i = 1:Nz
    for j = 1:Nr
        
        % set up the index
        n = inline('(i-1)*Nr + j','i','j','Nr');
        
        z_var(j,i) = z_initial + (i-1)*dz;
        r_var(j,i) = 0 + (j-1)*dr;
        
        z(i) = z_initial + (i-1)*dz;
        r(j) = 0 + (j-1)*dr;
        
        zeta = zeta_initial + (i-1)*dzeta;
        eta = 0 + (j-1)*deta;
        
        % set up the boundary conditions
        
        % Wall BC for R
        if(j == Nr  &&  (zeta_initial + (i-1)*dzeta) < 0)   % T0 at the wall
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            b_vec(n(i,j,Nr)) = 0;  % scaled variable - equivalent to T0
%         elseif(j == Nr && (zeta_initial + (i-1)*dzeta) < 1)             % T1 at the wall
%             A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
%             b_vec(n(i,j,Nr)) = zeta;  % scaled variable - equivalent to T1
        elseif(j == Nr)             % T1 at the wall
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            b_vec(n(i,j,Nr)) = 1; %zeta/zeta_final;  % scaled variable - equivalent to T1
        end
        % Centerline BC: dT/dr = 0
        if(j == 1)             % dT/dr = 0  --> polynomial interpolation
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            A_sparse(n(i,j,Nr),n(i,j+1,Nr)) = -4/3;
            A_sparse(n(i,j,Nr),n(i,j+2,Nr)) = 1/3;
            b_vec(n(i,j,Nr)) = 0;
        end
        
        % z-direction BCs
        if(i == 1  && j > 1 && j < Nr)
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            b_vec(n(i,j,Nr)) = 0;  % scaled variable - equivalent to T0
        elseif(i == Nz  && j > 1 && j < Nr)
            %A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            %b_vec(n(i,j,Nr)) = 1;  % scaled variable - equivalent to T1
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            A_sparse(n(i,j,Nr),n(i-1,j,Nr)) = -4/3;
            A_sparse(n(i,j,Nr),n(i-2,j,Nr)) = 1/3;
            b_vec(n(i,j,Nr)) = 0;
        end
        
        % now set up the equations for the internal region
        if(i > 1 && i < Nz  &&  j > 1 && j < Nr)
            b_vec(n(i,j,Nr)) = 0;
            
            A_sparse(n(i,j,Nr),n(i,j+1,Nr)) = eta/deta^2 + 1/(2*deta);
            A_sparse(n(i,j,Nr),n(i+1,j,Nr)) = eta/(Pe^2*dzeta^2);
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = -2*eta/deta^2 - 2*eta/(Pe^2*dzeta^2) - eta*(1-eta^2)/dzeta;
            A_sparse(n(i,j,Nr),n(i,j-1,Nr)) = eta/deta^2 - 1/(2*deta);
            A_sparse(n(i,j,Nr),n(i-1,j,Nr)) = eta/(Pe^2*dzeta^2) + eta*(1-eta^2)/dzeta;
            
        end
        
    end
end

%figure; spy(A_sparse)

theta_vec = A_sparse\b_vec;
%[theta_vec,flag,relres,iter] = gmres(A_sparse,b_vec,3,1e-12,10,[],[],[]);

resid = A_sparse*theta_vec - b_vec;

sum(resid)

theta = reshape(theta_vec,Nr,Nz);

Temp = 300 + theta*(100);

figure;
surf(z_var,r_var,theta); shading interp; colormap jet;

figure;
pcolor(z_var,r_var,Temp); colorbar; shading interp; colormap jet;

figure;contour(z_var,r_var,Temp,50); colorbar; shading interp; colormap jet;

figure;
plot(z,theta(0*(Nr-1)+1,:),'-',z,theta(0.2*(Nr-1)+1,:),'--',z,theta(0.4*(Nr-1)+1,:),':',...
    z,theta(0.6*(Nr-1)+1,:),'--',z,theta(0.8*(Nr-1)+1,:),':');