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

1 view (last 30 days)
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
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
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.

Pranav Verma on 15 Sep 2020

R2020a

### Community Treasure Hunt

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

Start Hunting!

Translated by