Spatial discretization has failed. Discretization supports only parabolic and elliptic equations.

1 view (last 30 days)
Mohammed Khan
Mohammed Khan on 10 Sep 2020
Answered: Pranav Verma on 15 Sep 2020
I am trying to solve two coupled Heat equations. One is a pde while the other is an ODE.
Is that why I am getting the error?:
Error using pdepe (line 293)
Spatial discretization has failed.
Discretization supports only parabolic and
elliptic equations, with flux term
involving spatial derivative.
Error in PBR1D>pdecalc (line 106)
sol = pdepe(m,@pdefun,@pdeic,@pdebc,z,t);
Error in PBR1D (line 11)
sol = pdecalc(h,P);
Equation information is found here:
CODE:
%% PBR cooling function
%% Contains both substrate and purge gas temp functions
%% Initial condition is from TMA half-reaction
%% NEGLECTS RADIAL HEAT TRANSFER
P = parameters();
P = bedPorosity(P);
[h,P] = heatTransferCoefficient(P);
sol = pdecalc(h,P);
function [P] = parameters()
P = zeros(1,40);
%% 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
P(4) = 889.2506; % [J/kgK] heat capacity of substrate
P(5) = 2330; %[kg/m3] density of silica
%
% P6 = E, P7 = Av
%
% skipped P(8) by accident
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 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
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
%% PDE Functions
function [sol] = pdecalc(h,P)
t = linspace(0,600,1200);
z = linspace(0,0.2,100); % length of reactor is 200mm
m = 0;
sol = pdepe(m,@pdefun,@pdeic,@pdebc,z,t);
%% Define: gam, alph, bet, C and A
function [c,f,s] = pdefun(z,t,T,dTdz) % PDE system to solve
% Parameter calculations
bet = (h*P(7))/(P(2)*P(3)*P(6));
alph = -P(1)/P(6);
gam = P(20);
A = P(4)*P(5)*(1-P(6));
C = h*P(7);
c = [1; 1];
f = [gam; 0] .*dTdz + [alph; 0];
s = [bet*(T(2)-T(1)); (-C/A)*(T(2)-T(1))];
end
function T0 = pdeic(z) % Initial conditions
T0 = [398.15; 548.2];
end
%% Define R
function [pl,ql,pr,qr] = pdebc(zl,Tl,zr,Tr,t) % Boundary conditions
R = P(21);
T0 = P(22);
%pl = [(-1/R)*Tl(1) + (1/R)*T0; 0];
pl = [Tl(1)-T0; 0];
ql = [-R; 0];
pr = [0; 0];
qr = [1; 0];
end
end
  1 Comment
Mohammed Khan
Mohammed Khan on 10 Sep 2020
I have debugged the issue and found that the issue is due to the boundary conditions I defined. I would really appreciate any help or tips in the proper formulation/implementation.

Sign in to comment.

Answers (1)

Pranav Verma
Pranav Verma on 15 Sep 2020

Community Treasure Hunt

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

Start Hunting!

Translated by