% 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 p1_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 = 5*4 + 1;
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);


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;
        
        v_z = 2*U*(1-(r(j)/R)^2);

        
        % set up the boundary conditions
        
        % Wall BC for R
        if(j == Nr  &&  (z_initial + (i-1)*dz) < 0)   % T0 at the wall
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            b_vec(n(i,j,Nr)) = 300;  % scaled variable - equivalent to T0
        elseif(j == Nr)             % T1 at the wall
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = 1;
            b_vec(n(i,j,Nr)) = 400; %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)) = 300;  % 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)) = 1/dr^2 + 1/(2*r(j)*dr);
            A_sparse(n(i,j,Nr),n(i+1,j,Nr)) = 1/(dz^2);
            A_sparse(n(i,j,Nr),n(i,j,Nr)) = -2/dr^2 - 2/(dz^2) - v_z/alpha/dz;
            A_sparse(n(i,j,Nr),n(i,j-1,Nr)) = 1/dr^2 - 1/(2*r(j)*dr);
            A_sparse(n(i,j,Nr),n(i-1,j,Nr)) = 1/(dz^2) + v_z/alpha/dz;
            
        end
        
    end
end

Temp_vec = A_sparse\b_vec;
%[Temp_vec,flag,relres,iter] = gmres(A_sparse,b_vec,5,[],100,[],[],[]);

resid = A_sparse*Temp_vec - b_vec;

sum(resid)

Temp = reshape(Temp_vec,Nr,Nz);

figure;
surf(z_var,r_var,Temp); shading interp; colormap jet;
axis([z_initial z_final 0 R 300 400]);
xlabel('Length (meters)'); ylabel('Radius (cm)'); zlabel('Temperature (K)');

figure;
pcolor(z_var,r_var*100,Temp); colorbar; shading interp; colormap jet;
axis([z_initial z_final 0 R*100]);
xlabel('Length (meters)'); ylabel('Radius (cm)'); title('Temperature Field (K)');

figure;
[C,h] = contour(z_var,r_var*100,Temp,10); colorbar; shading interp; colormap jet;
axis([z_initial z_final 0 R*100]);
%set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2)
xlabel('Length (meters)'); ylabel('Radius (cm)'); title('Temperature Field (K)');

figure;
plot(z,Temp(0*(Nr-1)+1,:),'-',z,Temp(0.2*(Nr-1)+1,:),'--',z,Temp(0.4*(Nr-1)+1,:),':',...
    z,Temp(0.6*(Nr-1)+1,:),'--',z,Temp(0.8*(Nr-1)+1,:),':');
axis([z_initial z_final 300 400]);  
xlabel('Length (meters)'); ylabel('Temperature (K)'); title('Temperature at various r-values');
legend('r = 0 cm','r = 1 cm','r = 2 cm','r = 3 cm','r = 4 cm');