Need help using pdepe

Equation(1)

This is the PDE used to calculate the mole fraction of different species (Y) with input gas flow rate (F) and growth temperature (T). The code provided is to for Y of one species which is H2.

Range of F is 1.667e-7 to 1.667e-6

Range of T is 773 to 1473

The output that I need is

1. Graph of Y vs T at constant F
2. Graph of Y vs F at constant T

My problem is that I couldn’t solve the given PDE using MATLAB.

In order to use pdepe we have to obtain the coefficient c, f, s and m:

My Equation(1) has one extra term compared the one presented in the documentation

Here is my initial condition:

Y (0,0) = 0.2, Y_H2 = 0.2

And here is my boundary condition:

My system is a slab with x boundaries to describe the width of the slab

x1 = 0 and x2 = 0.02

x1 and x2 are the limits of the slab which is the length of the slab

Below is my current code please help if you can

%% setting the values of m = the molecular weight of the species(kg/kmol)
H2_m = 2.016;
%% setting the values of rho = the Density of the species(kg/m3)
H2_rho = 0.08988;
%% setting the values of D = the species mass diffusion coefficient(m2/s)
H2_D = 1.233464888;
%% Area of the reactor
L = 0.02; % Length in meter
W = 0.15; % Width in meter
Area = L*W; % Area of the reactor = (Length x Width) m2
A = [1.00000000000000e+18	9.20000000000000e+16	22000	90000000000000.0	690000000000000	1.00000000000000e+18	70000000000000.0	30000000000000.0	150000000000000	540	125000000000000	3.36000000000000e-07	18000000000000.0	40000000000000.0	2.54500000000000e+19	90000000000.0000	90000000000.0000	1100000000000.00	28000000	28000000	28000000	3.20000000000000e+30	3.20000000000000e+30	3.20000000000000e+30];
B = [-1	-0.600000000000000	3	0	0	-1.56000000000000	0	0	0	3.50000000000000	0	6	0	0	0	2	2	0	2	2	2	0	0	0];
Ea = [0	0	8750	15100	82469	0	0	0	0	5210	8000	1692	76000	0	19379	5000	5000	7300	7700	7700	7700	71945	71945	71945];
vik = [-1	-1	-1	-1	-1	-1	1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-3	-1	-2];
R = 1.987;
n = length(A);
F = linspace((1.667*10^-7),(1.667*10^-6),n); % Inlet flow rate in (m3/s)
T = linspace(773,1473,n); % Growth Temperature in K
R_reaction_array = zeros(1,n);

%% Start Caluclations
for i = 1:n
R_reaction_array(i) = vik(i)*(  A(i)*(T(i)^B(i))*exp(-Ea(i)/(R*T(i)))   );
end
%% ________________________________________
function [c,f,s] = pdefun(x,t,Y,dYdx)
c = 1+F/Area;
D = H2_D;
f = -D*(dYdx);
s = (H2_m/H2_rho)*sum(R_reaction_array);
end
%% ________________________________________
function Y0 = icfun(x)
Y0 = 0.2;
end
%% ________________________________________
function [pL,qL,pR,qR] = bcfun(xL,YL,xR,YR,t)
pL = YL;
qL = 0;
pR = YR - 1;
qR = 0;
end
%% ________________________________________

x = linspace(0,0.02,25);
t = linspace(0,10,25);
m = 0;
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);

[Torsten]

function [c,f,s] = pdefun(x,t,Y,dYdx)
  c = 1.0;  
  f = H2_D*dYdx;
  s = -F/Area*dYdx + (H2_m/H2_rho)*sum(R_reaction_array);
end

[Abdullah Alhebshi]

@Torsten I get these warnings

Warning: Failure at t=0.000000e+00.  Unable to meet integration tolerances without reducing the step size below
the smallest value allowed (7.905050e-323) at time t.

Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.

[Torsten]

Start simple.

Does the following code work ?

function main
%% setting the values of m = the molecular weight of the species(kg/kmol)
H2_m = 2.016;
%% setting the values of rho = the Density of the species(kg/m3)
H2_rho = 0.08988;
%% setting the values of D = the species mass diffusion coefficient(m2/s)
H2_D = 1.233464888;
%% Area of the reactor
L = 0.02; % Length in meter
W = 0.15; % Width in meter
Area = L*W; % Area of the reactor = (Length x Width) m2
A = [1.00000000000000e+18	9.20000000000000e+16	22000	90000000000000.0	690000000000000	1.00000000000000e+18	70000000000000.0	30000000000000.0	150000000000000	540	125000000000000	3.36000000000000e-07	18000000000000.0	40000000000000.0	2.54500000000000e+19	90000000000.0000	90000000000.0000	1100000000000.00	28000000	28000000	28000000	3.20000000000000e+30	3.20000000000000e+30	3.20000000000000e+30];
B = [-1	-0.600000000000000	3	0	0	-1.56000000000000	0	0	0	3.50000000000000	0	6	0	0	0	2	2	0	2	2	2	0	0	0];
Ea = [0	0	8750	15100	82469	0	0	0	0	5210	8000	1692	76000	0	19379	5000	5000	7300	7700	7700	7700	71945	71945	71945];
vik = [-1	-1	-1	-1	-1	-1	1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-1	-3	-1	-2];
R = 1.987;
n = length(A);
F = linspace((1.667*10^-7),(1.667*10^-6),n); % Inlet flow rate in (m3/s)
T = linspace(773,1473,n); % Growth Temperature in K
R_reaction_array = zeros(1,n);
%% Start Caluclations
for i = 1:n
    R_reaction_array(i) = vik(i)*(  A(i)*(T(i)^B(i))*exp(-Ea(i)/(R*T(i)))   );
end
x = linspace(0,0.02,25);
t = linspace(0,10,25);
m = 0;
sol = pdepe(m,@(x,t,Y,dYdx)pdefun(x,t,Y,dYdx,H2_D,F(1),Area),@icfun,@bcfun,x,t);
u = sol(:,:,1);
surf(x,t,u)
xlabel('x')
ylabel('t')
zlabel('u(x,t)')
view([150 25])
end
%%____________________________________
function [c,f,s] = pdefun(x,t,Y,dYdx,H2_D,F,Area)
c = 1.0;
f = H2_D*dYdx;
s = -F/Area*dYdx;
end
%% ________________________________________
function Y0 = icfun(x)
Y0 = 0.2;
end
%% ________________________________________
function [pL,qL,pR,qR] = bcfun(xL,YL,xR,YR,t)
pL = YL;
qL = 0;
pR = YR - 1;
qR = 0;
end
%% ________________________________________

[Abdullah Alhebshi]

@Torsten Yes it works.

But s used to be

s = -F/Area*dYdx + (H2_m/H2_rho)*sum(R_reaction_array);

you changed it to

s = -F/Area*dYdx

by removing this term

(H2_m/H2_rho)*sum(R_reaction_array)

May I know why it wroked after you removed it ?

and is it correct to remove it? Because its a whole term in the PDE in Equation(1)

[Torsten]

But why do you have to add all the reaction terms which are specified at different temperatures ?

Are you sure this is correct ?

[Abdullah Alhebshi]

@Torsten Because this PDE describes a process in which multiple reactions are taking place. Thus, the effect of temperature must be accounted for, in order to study the effect of reaction temperature on the final mole fraction (Y).

it is necessary to account for the effect of temperature

[Torsten]

But different temperatures for one study case?

[Abdullah Alhebshi]

@Torsten Yes