Error using pdepe: Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative

I get this error trying to solve a system of PDEs, and I do not know if such system is solvable with 'pdepe'. The equations are:

energy balance equation, and:

mass balance equations. Initial conditions:

Boundary conditions:

The variables in which I am solving are however T,, knowing that:

, I can write initial condition for , and I have also a value for .

I tried writing the equations in the form required by 'pdepe'. Also, note that many parameters depends on the vector , IN particular, the term 'v' includes a partial spatial derivative w.r.t. P, which i transformed in a spatial derivative w.r.t. u(1) and u(2) with the above formula. Note that many constants are imported with the file.m including the initial conditions that are: This is the code of the 3 functions used inside the pde call (sorry is a bit long and messy):

function [c,f,s] = pdefun(x,t,u,dudx)
        % load data
        run('HM_data.m');

        % compute epsilon , F and dfdt
        F = (hm.rho_s_in-u(3))/(u(3)*hm.tau_p - hm.rho_s_in*hm.w_max);
        epsilon = 1 - (1 - hm.epsilon_0)*((1 + hm.tau_p*F)/(1 + hm.tau_a*F));
        P = u(2)*u(1)*hm.R_gas/hm.MH2;
        P_eq_a = (hm.C0_a*hm.C1_a*(F)^(hm.C2_a)/(1+hm.C1_a*(F)^(hm.C2_a)) + hm.C3_a*F + exp(hm.C4_a*(F-hm.C5_a)))*exp(-hm.K_a*((1/u(1))-(1/303))); 
        k_a = (hm.kappa_a/(1 - epsilon)) * exp(-hm.E_a/(hm.R_gas*u(1))) * log(P/P_eq_a);
        dfdt = k_a*(1-F);
       

        % compute diagonal matrix c
        Cps = 6000*(3.1*hm.R + 10.04*hm.x_max*F)/(hm.Ms + 6*hm.x_max*F);
        rho_C_eff =  epsilon*u(2)*hm.Cpg + (1- epsilon)*u(3)*Cps;
        c = [rho_C_eff; epsilon; 1-epsilon];

        % COMPUTE VECTOR f
        % velocity
        Dp = hm.D_in*(1+ F*hm.tau_p)^(1/3);
        Kp = (Dp^2)*(epsilon^3)/(150*(1-epsilon)^2);
        v_cost = -Kp*hm.R_gas/(hm.mu_g*hm.MH2);
        v = v_cost*(dudx(1)*u(2) + dudx(2)*u(1));

        % effective thermal conductivity
        N = 3.08/epsilon - 1.13;
        Fn = 4/3 * hm.E_prime * hm.R^(1/2) * hm.dv^(1.5);
        Hv = hm.c1*hm.dv^(hm.c2);
        Rs = 0.565*Hv*hm.dv/(hm.ks*Fn);
        a_H = (0.75*Fn*hm.R/hm.E_prime)^(1/3);
        a_LH = [];
        if epsilon <= 0.47 && epsilon >= 0.01
                a_LH = 1.605/sqrt(epsilon);
        elseif epsilon > 0.47 && epsilon <= 1
                a_LH = 3.51 - 2.51*epsilon;
        else 
                error('error epsilon');
        end
        a_L = a_LH*a_H;
        R_L = 1/(2*hm.ks*a_L);

        Pmax = 2/pi*hm.E_prime*(hm.dv/hm.R)^(0.5);
        Hc = hm.c1*(1.62*hm.dv)^(hm.c2);
        a1 = erfcinv(2*Pmax/Hc);
        a2 = erfcinv(0.03*Pmax/Hc) - a1;
        coef_a = 2*hm.b/Dp;
        coef_c = -(6*(hm.gamma-1)/(9*hm.gamma-5))* (hm.kg_ref*hm.MH2/(u(1)*u(2)*hm.R_gas))*(hm.MH2*u(1)/(2*hm.kb))^(0.5);
        
        l_m = (-1 + sqrt(1 - 4*coef_a*coef_c))/(2*coef_a);
        kg = hm.kg_ref/(1+2*hm.b*l_m/Dp);
        M = (((2-hm.alpha_T1)/hm.alpha_T1)+((2-hm.alpha_T2)/hm.alpha_T2))*(2*hm.gamma/(1+hm.gamma))*(1/hm.Pr)*l_m;
        R_g = sqrt(2)*hm.sigma*a2/(pi * kg * a_L^2 * log(1+ a2/(a1+M/(sqrt(2)*hm.sigma))));
        L = (hm.gamma+1)*3*Dp/((9*hm.gamma-5)*4*l_m*sqrt(pi));
        R_G = 1/(2*pi*kg*Dp*(0.5*log(1+L) + log(1+sqrt(L)) + 1/(1+sqrt(L)) - 1));
        R_mic = Rs*R_g/(Rs+R_g);
        R_c_inv = 1/(R_mic+R_L) + 1/R_G;

        k_eff = N*(1-epsilon)*R_c_inv/(pi*Dp);

        f = [k_eff*dudx(1); -u(2)*v; 0];
         
        % compute vector s
        m_dot_a = (1-epsilon)*(hm.rho_sat-u(3))*dfdt;
        s = [u(2)*hm.Cpg*v*dudx(1) + m_dot_a*hm.delta_H ; -m_dot_a+hm.phi_abs; m_dot_a];

end


function u0 = icfun(x)
        run('HM_data.m');
        u0 = hm.ic;
end

function [pl,ql,pr,qr] = bcfun(xl, ul, xr, ur, t)
        run('HM_data.m');
        Nu_abs = 0.3+((0.62*(hm.Re_a^(0.5))*(hm.Pr_a^(1/3)))/(1+(0.4/hm.Pr_a)^(2/3))^(1/4)*((1+(hm.Re_a/282000)^(5/8))^(4/5)));
        h_f = Nu_abs*hm.kf/hm.D_tank;
        pl = [0 ; 0; 0];
        ql = [1; 0; 0];
        pr = [h_f*(ur(1)-hm.Tf_a); 0; 0];
        qr = [-1; 0; 0];
end

Answer 1

[Torsten]

"pdepe" is a solver for parabolic-elliptic equations. Thus all equations should have a second-order derivative term in space and boundary conditions on the left and on the right.

This is true for your energy equation, but not for your mass balance equations.

You set

pl = [0 ; 0; 0];
ql = [1; 0; 0];
pr = [h_f*(ur(1)-hm.Tf_a); 0; 0];
qr = [-1; 0; 0];

thus no condition at either end of the interval for rho_g and rho_s. This results in 0 = 0 for the solver and it exits because no information is given.

With some tricks, you could set artificial conditions, but the results will not be trustworthy.

Summarizing: "pdepe" is not suited for your problem.

[Lorenz.zz]

Thank you very much for your response. Could you please give me some information on these 'tricks', or otherwise some other methods to use? Are these equations even solvable in Matlab? Also, I imagine that a method which solves all of the equations would be needed since they are coupled.

In the study work that I am following, they transformed the equations in this way in order to solve them, but I do not understand how:

for the two mass balance equations, saying it is solved with the method of characteristics, and :

for the energy equation, where , and new boundary conditions:

solved by combining separation of variables and superposition , summing homogeneous and steady-state solutions. But I do not understand which are the terms , and I could I solve this in Matlab. Thank you.

[Torsten]

The first thing is that you work in cylindrical coordinates - thus you have to set m = 1 in the call to "pdepe". Maybe you did this - you didn't include your call to the solver.

The equation for rho_g needs an inflow boundary condition at one end of the domain. You didn't set that in your code. Once you've set this condition, I wonder where the mass enters and leaves the domain since you defined the domain from r = 0 to r = R.

The "trick" for all equations where no boundary condition is known is to set p = 0, q = 1 in the boundary condition part.

Nevertheless, since the results won't be trustworthy in my opinion, I suggest discretizing the spatial derivatives while leaving the temporal derivatves in their continuous form and solve the resulting system of ordinary differential equations using ode15s. Look up "method-of-lines" for more details.

[Lorenz.zz]

Thank you for your response.

Yes I have set m=1 and called the function:

t_max = 3600;
t = linspace(0, t_max, 10);
x = linspace(0, hm.R, 10);


% computation of the solution of the PDEs system
m = 1;
sol = pdepe(m, @pdefun, @icfun, @bcfun, x, t);

% Extract solutions
T = sol(:,:,1);
rho_g = sol(:,:,2);
rho_s = sol(:,:,3);

The boundary conditions of the densities are not specified, by logic for cylindrical symmetry the gradients density should be 0 at r=0 . At the border r=R, since is the density of the gas that can enter the tank, the condition should be the entering gas flux (which by the way is the therm in the equation)? I am not sure since the inlet site is not at r=R but at one end of the cylinder. At r=R the only thing present is the contact of the tank with the external circulating fluid, but it does not effect masses of the solid nor gas, only temperature. While change indirectly with the entrance of the gas since there is a process of absorption within the solid, so i wouldn't know. My question is, is it possible that the boundary conditions for the mass balance equations were not given since they are considered '0' at both ends? I am trying to implement 'ode15s' with null boundary conditons but i do not know if it is the right path.

[Torsten]

I am not sure since the inlet site is not at r=R but at one end of the cylinder.

You have a radial velocity v in your equation for rho_g, not an axial. In case of axial flow, you have a 2d-problem in space with independent variables z and r.