Packed Bed heat transfer UPDATED

%% BEP MST Mohammed
%% PBR cooling function Method of lines
%% Contains both substrate and purge gas temp functions
%% Initial condition is from TMA half-reaction
%% NEGLECTS RADIAL HEAT TRANSFER
%% Dealing with PDE and ODE coupling

P = parameters();
P = bedPorosity(P);
[h,P] = heatTransferCoefficient(P);

% Defining Initial conditions
% Tf initial = 398.15K & Ts initial = 548.3K
N=10;
Tf(N+1,1)=0;
Tf(:,:)=398.15;
Ts(N+1,1)=0;
Ts(:,:) = 548.3;
y = zeros(N*2,1);
y(1:N)=Tf(1:N);
y(N+1:2*N)=Ts(1:N);

y0=y;

tSpan = linspace(0,60,120000);

tic
[t,y] = ode15s(@(t,y) heatPBR(t,y,h,P), tSpan, y0);
toc
% plot results
plot(t,y);


function [P] = parameters()
P = zeros(1,100);
%% All important variables are stored in the P array

% all important parameters for fluid PDE
P(1) = 1.0;         % [m/s] Velocity of N2
P(2) = 0.85;        % [kg/m3] Density of N2 | 400K,1bar
P(3) = 1044;        % [J/kgK] Heat capacity N2 | 400K

% all important parameters for Solid PDE
P(4) = 889.2506;    % [J/kgK] heat capacity of substrate
P(5) = 2330;        %[kg/m3] density of silica
%
% P6 = E, P7 = Av
%
P(8) = 1.0;                     % [m/s] velocity of N2
P(9) = 3E-3;                    % [m] diameter of silica particle
P(10) = 0.01;                   % [kg] Mass of particle
P(11) = 0.25*pi*0.2*(0.025^2);  % [m3] reactor volume
                                % Ruud: Length = 200mm, diam = 25mm.
P(12) = pi*(P(9)^2);            % [m2] Surface area of one particle
% P13 = Volume of all particles
P(14) = 2.12E-5;                % [Pa*s] Viscosity of N2
P(15) = 37E-6;                  % [m2/s] thermal diffusivity N2 at 1 Bar!!
P(16) = 32.81E-3;               % [W/mK] thermal conductivity N2 at 400K
%
% P17 = Nu, P18 = Re, P19 = Pr
%
P(20) = P(16)/(P(2)*P(3));      % [] Gamma
P(21) = P(16)/(P(8)*P(2)*P(3)); % R variable in BC's
P(22) = 398.15;                 % [K] initial value N2
P(23) = 1.4;                    % [W/mK] thermal conductivity silica
end

% calculate porosity of bed
function [P] = bedPorosity(P)
d = P(9);                   % [m] diameter of silica particle
rho = P(5);                 % [kg/m3] density of silica particle
m = P(10);                  % [kg] mass of particles
vs = m/rho;                 % [m3] volume of particles
v = P(11);                  % [m3] Volume reactor
                            % Ruud: Length = 200mm, diam = 25mm.
vVoid = v-vs;
E = vVoid/v;                % porosity

Ap = P(12);                             % [m2] Surface area of one particle
Vsp = 0.25*(pi*(0.025^2))*0.2*(1-E);    % [m3] Volume of all particles;                            % [m3] volume all substrate
Av =  Ap/Vsp;                           % [m2/m3] this Av is calculated as:
                                        % Area 1 sphere / Volume all spheres
                                        % Area 1 sphere: pi*(3E-3)^2
                                        % Volume all spheres: 0.25*(pi*(0.025^2))*0.2*(1-E)
P(6) = E;
P(7) = Av;
P(13) = Vsp;
end

function [h,P] = heatTransferCoefficient(P)
rho = P(2);         % [kg/m3] density of N2
v = P(1);           % [m/s] velocity ACTUALLY related to phiV!!!
d = P(9);           % [m] diameter of 1 silica particle
mu = P(14);         % [Pa*s] Viscosity of N2
a = P(15);          % [m2/s] thermal diffusivity N2 at 1 Bar!!
nu = mu/rho;        % [m2/s] kinematic viscosity
lambda = P(16);     % [W/mK] thermal conductivity at 400K

Re = (rho*v*d)/mu;      % Reynolds number. Depends on density, velocity,
                        % viscosity of fluid. d is diameter of substrate
Pr = nu/a;              % Prandtl number, nu is kinematic viscosity of
                        % fluid, a is thermal diffusivity
E = P(6);

%% Gnielinski method
% E < 0.935, Pr>0.7, Re/E<2E4
Nut = (0.037*((Re/E)^0.8)*Pr)/(1+2.433*((Re/E)^-0.1)*(Pr^(2/3)-1));
Nul = 0.664*(Pr^(1/3))*((Re/E)^0.5);
Nusp = 2+((Nut+Nul)^0.5);
fE = (1+1.5*(1-E));
Nu = Nusp*fE;           % Nusselt number | not to be confused with 'nu'

h = (Nu*lambda)/d;      % Heat transfer coefficient
P(17) = Nu;
P(18) = Re;
P(19) = Pr;
end

%% Method of Lines
function dTdt = heatPBR(t,y,h,P)
% Parameters
alpha= -P(1)/P(6);
gamma= P(20);
beta= (h*P(7))/(P(2)*P(3)*P(6));
omega= P(23)/(P(2)*P(3)*(1-P(6)));

% Getting temperatures
N = length(y);
M = length(y)/2;
dz = 0.2/N;                 % length of reactor is 200mm
dz2 = dz^2;
Tf = zeros(M,1);
Ts = zeros(M,1);
Tf(1:M)=y(1:M);
Ts(1:M)=y(M+1:2*M);

% Define dT/dt
dTfdt=zeros(M,1);
dTsdt=zeros(M,1);

% Boundary conditions
R = P(21);
T0 = P(22);
a = P(1)*P(2)*P(3);

for i=1:M
    if i==1
        dTfdt(i) = alpha*((-a/P(16))*(T0-Tf(i))) + gamma*((-2*Tf(i)-2*a*dz*(T0-Tf(i))/dz2)) + beta*(Ts(i)-Tf(i));
    elseif i==M
        dTfdt(i) = gamma*((2*Tf(i-1)-2*Tf(i))/dz2) + beta*(Ts(i)-Tf(i));
    else
        dTfdt(i) = (alpha/(2*dz))*(Tf(i+1)-Tf(i-1)) + (gamma/dz2)*(Tf(i+1)-2*Tf(i)+Tf(i-1)) + beta*(Ts(i)-Tf(i));
    end
end

for i=1:M
    dTsdt(i) = omega*(Ts(i)-Tf(i));
end

dTdt=zeros(2*M,2);
dTdt=[dTfdt;dTsdt];
end