%2D Finite-difference transient conduction calculator
%Greg Teichert & Kyle Halgren
%This program calculates the transient conduction of a rectangular cross
%section.
%
%
%                     h(2), T_inf(2)
%                   _________________
%                  |                 |
%   h(1), T_inf(1) |                 | a   h(3), T_inf(3)
%                  |_________________|
%                          b
%                     h(4), T_inf(4)

%%Inputs
%Convective heat transfer coefficients
%For a side to be completely insulated, set h = 0.
%For a side to be a constant temperature, set h to be very high
%The program assumes that the convective coefficients of the corners are
%the average of the two adjecent sides
h = zeros(4,1);%this is to set up thevector of h values. DO NO CHANGE.

%Userdefined h values
h(1) = 10;
h(2) = .1;
h(3) = 10;
h(4) = .1;

%Boundary conditions
%Userdefined T infinity values in kelvin
T_inf(1) = 293;
T_inf(2) = 293;
T_inf(3) = 353;
T_inf(4) = 353;

%Initial condition (assume uniform initial temperature)
%Userdefined initial temperature value
T_0 = 573;

%Material properties
%Userdefined material values
k = .08;
rho = 7480;
c_p = .460;

%Userdefined physical variables
a = 1; %height of cross section
b = 1.3; %width of cross section
t = 3600; %time at which results are given

%%Calculations
%alpha (thermal diffusivity)
alpha = k/(rho*c_p);


%delta_x (space increment, equal to delta_y)
delta_x = gcd(a*100000,b*100000)/100000;
while delta_x > min(a,b)/20
    delta_x = delta_x/2;
end

%Biot numbers
Bi = zeros(4,1);
Bi = h*delta_x/k;

%delta_t (time increment)
Bi_max = max(Bi);
delta_t1 = 0.25*delta_x^2/alpha;
delta_t2 = delta_x^2/(2*alpha*(2 + Bi_max));
delta_t3 = delta_x^2/(4*alpha*(1 + Bi_max));
delta_t = min([delta_t1 delta_t2 delta_t3]);

%Fourier number
Fo = alpha*delta_t/delta_x^2;

%%Create matrices for solution
%Creation of two constants that are frequently used in the code; they are
%the number of nodes in a column and row of nodes, respectively.
a_step = a/delta_x + 1;
b_step = b/delta_x + 1;
totalnodes = a_step*b_step;

%Master temperature matrix: Each column contains the nodal temperatures of
%the cross section for one time step. The columns are made made going down
%the column of nodes in the cross section, from left to right.
T = zeros(totalnodes,round(t/delta_t) + 1);

%Initial condition
T(1:totalnodes) = T_0*ones(totalnodes,1);



for k = 2:(round(t/delta_t) + 1)
    %Setting up the coefficient matrix and constant matrix
    coeff = zeros(totalnodes,totalnodes);
    const = zeros(totalnodes,1);
    for i = 1:totalnodes
        %Top left corner
        if i == 1
            coeff(i,i + a_step) = 2*Fo;
            coeff(i,i + 1) = 2*Fo;
            const(i,1) = 4*Fo*(Bi(1)*T_inf(1) + Bi(2)*T_inf(2))/2;
            coeff(i,i) = 1 - 4*Fo - 4*Fo*(Bi(1) + Bi(2))/2;
        %Top right corner
        elseif i == (totalnodes - a/delta_x)
            coeff(i,i - a_step) = 2*Fo;
            coeff(i,i + 1) = 2*Fo;
            const(i,1) = 4*Fo*(Bi(2)*T_inf(2) + Bi(3)*T_inf(3))/2;
            coeff(i,i) = 1 - 4*Fo - 4*Fo*(Bi(2) + Bi(3))/2;
        %Bottom left corner
        elseif i == a_step
            coeff(i,i + a_step) = 2*Fo;
            coeff(i,i - 1) = 2*Fo;
            const(i,1) = 4*Fo*(Bi(4)*T_inf(4) + Bi(1)*T_inf(1))/2;
            coeff(i,i) = 1 - 4*Fo - 4*Fo*(Bi(4) + Bi(1))/2;
        %Bottom right corner
        elseif i == totalnodes
            coeff(i,i - a_step) = 2*Fo;
            coeff(i,i - 1) = 2*Fo;
            const(i,1) = 4*Fo*(Bi(3)*T_inf(3) + Bi(4)*T_inf(4))/2;
            coeff(i,i) = 1 - 4*Fo - 4*Fo*(Bi(3) + Bi(4))/2;
        %Left side
        elseif 1 < i  && i < a_step
            coeff(i,i + a_step) = 2*Fo;
            coeff(i,i - 1) = Fo;
            coeff(i,i + 1) = Fo;
            const(i,1) = 2*Fo*Bi(1)*T_inf(1);
            coeff(i,i) = 1 - 4*Fo - 2*Bi(1)*Fo;
        %Right side
        elseif (totalnodes - a/delta_x) < i && i < totalnodes
            coeff(i,i - a_step) = 2*Fo;
            coeff(i,i - 1) = Fo;
            coeff(i,i + 1) = Fo;
            const(i,1) = 2*Fo*Bi(3)*T_inf(3);
            coeff(i,i) = 1 - 4*Fo - 2*Bi(3)*Fo;    
        %Bottom
        elseif rem(i/a_step,2) == 0 || rem(i/a_step+1,2) == 0
            coeff(i,i - 1) = 2*Fo;
            coeff(i,i + a_step) = Fo;
            coeff(i,i - a_step) = Fo;
            const(i,1) = 2*Fo*Bi(4)*T_inf(4);
            coeff(i,i) = 1 - 4*Fo - 2*Bi(4)*Fo;    
        %Top
        elseif rem((i - 1)/a_step,2) == 0 || rem((i - 1)/a_step+1,2) == 0
            coeff(i,i + 1) = 2*Fo;
            coeff(i,i + a_step) = Fo;
            coeff(i,i - a_step) = Fo;
            const(i,1) = 2*Fo*Bi(2)*T_inf(2);
            coeff(i,i) = 1 - 4*Fo - 2*Bi(2)*Fo;    
        %Inside
        else
            coeff(i,i + a_step) = Fo;
            coeff(i,i - a_step) = Fo;
            coeff(i,i + 1) = Fo;
            coeff(i,i - 1) = Fo;
            coeff(i,i) = 1 - 4*Fo;
        end
    end
    T(:,k) = const + coeff*T(:,k-1);
end

%%Create the final temperature distribution matrix of the cross section

T_distribution  = zeros(a_step,b_step);

for j = 1:b_step
    T_distribution(:,j) = T(((j - 1)*a_step + 1):j*a_step,round(t/delta_t) + 1);
end

%Plot T_distribution
x = 0:delta_x:b;
y = 0:delta_x:a;

figure(1)
clf
surf(y,x,T_distribution')


xlabel('a (m)')
ylabel('b (m)')
zlabel('Temperature (K)')
title('Transient conduction (the origin of the plot is the top left corner of the cross section)')