%'HT_1D_planar.m'. 
%One-dimensional heat transfer in planar geometry by finite differences. 
%Equation to solve: rho*c*DT/Dt=k*D2T/Dx2+phi-rho*c*v*DT/Dx
clear;close all;%Default data applicable to Problem P-11.53
%Data:          %Change data and see corresponding results
D=0.000;        %Diametre or characteristic transversal dimension [m]
L=0.100;        %Length along heat path [m]
A=0.15*1e-3;    %Area. Direct input, or computed from diameter as A=pi*D^2/4 [m2]
p=2*0.15;       %Perimeter. Direct input, or computed from diameter as p=pi*D [m]
k=0.25;         %Conductivity [W/(m·K)]
c=700;          %Thermal capacity [J/(kg·K)]
rho=1000;       %Density [kg/m3]
h=2;            %Convective coefficient [W/(m^2·K)], transversal
eps=0;          %Emissivity, transversal
sigma=5.67e-8;  %S-B constant
N=20;           %# of elements along L
M=400;          %# of time steps
Tinf=288;       %Ambient temperature [K]
tsim=10000;     %Total simulation time [s]
Troot=303;      %Root temperature [K] in case it is fixed
Qroot=0;        %Root heat transfer power [W] in case it is fixed
phi=5/(A*L);    %Volumetric dissipation [W/m^3] if any

%Inicialization:
a=k/(rho*c);            %Diffusivity [m^2/s]
Dx=L/N;                 %Element width
X=linspace(0,L,N+1);    %Node position list (equispaced)
Dt=tsim/M;              %Time step (you might fix it instead of tsim)
t=linspace(0,tsim,M)';  %Time vector
Fo=a*Dt/(Dx*Dx)         %Fourier's number
Bi=h*p*Dx/(k*A/Dx)      %Biot's number
%Check for stability of the explicit finite difference method
disp(['Stability requires 1-Fo*(2+Bi)<0. It actually is =',num2str(1-Fo*(2+Bi))])
if 1-Fo*(2+Bi)<0 disp('This is unstable; increase number of time steps'), end
T=Tinf*ones(M,N+1); 	%Temperature-matrix (times from 1 to M, and positions from 1 to N+1)
 
%Iteration:   
for it=2:M
    for i=2:N
        T(it,i)=T(it-1,i)+Fo*(T(it-1,i+1)-2*T(it-1,i)+T(it-1,i-1))+Fo*Bi*(Tinf-T(it-1,i))+phi*Dt/(rho*c);
    end
    %Boundory condition in node 0 (one of the following two):
    T(it,1)=Troot;      %if Troot is fixed
    %T(it,1)=T(it,2)-Qroot*Dx/(2*k*A);   %if Qroot is fixed
    %Boundory condition in node N (one of the following two):
    T(it,N+1)=Troot;    %if Troot is fixed
    %T(it,N+1)=T(it,N)-Qroot*Dx/(2*k*A);   %if Qroot is fixed
end
 
%Presentation of results
subplot(2,1,1);plot(t,T(:,1:N/10:N+1));xlabel('t [s]'),ylabel('T [K]');title('T(t,x) vs. t at several locations')
subplot(2,1,2);plot(X,T(1:M/100:M,:));xlabel('X [m]'),ylabel('T [K]');title('T(t,x) vs. X at several times')