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

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

[Mohammed Khan]

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.

Answer 1

[Pranav Verma]

Hi Mohammed,

Please refer to the below posts on the similar lines.

https://www.mathworks.com/matlabcentral/answers/290897-is-matlab-able-to-solve-my-system-of-coupled-ode-and-pde-equations

https://www.mathworks.com/matlabcentral/answers/168832-how-to-solve-pde-and-ode-together

https://www.mathworks.com/matlabcentral/answers/510000-how-do-i-solve-a-coupled-system-of-pde-and-ode

Thanks