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

20 views (last 30 days)
Lorenz.zz on 10 Jul 2024 at 15:00
Edited: Torsten on 12 Jul 2024 at 9:30
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)
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

Torsten on 10 Jul 2024 at 15:23
Edited: Torsten on 10 Jul 2024 at 15:24
"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 on 12 Jul 2024 at 7:57
Edited: Lorenz.zz on 12 Jul 2024 at 7:59
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 on 12 Jul 2024 at 9:30
Edited: Torsten on 12 Jul 2024 at 9:30
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.