Need help using pdepe

3 views (last 30 days)
Abdullah Alhebshi
Abdullah Alhebshi on 22 Dec 2021
Commented: Abdullah Alhebshi on 29 Dec 2021
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);
  8 Comments

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by