Error using pdepe: Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative
28 views (last 30 days)
Show older comments
Lorenz.zz
on 10 Jul 2024 at 15:00
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:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731738/image.jpeg)
energy balance equation, and:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731743/image.jpeg)
mass balance equations. Initial conditions:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731748/image.jpeg)
Boundary conditions:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731753/image.jpeg)
The variables in which I am solving are however T,
, knowing that:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731758/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731763/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731768/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731773/image.png)
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):
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731778/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1731783/image.png)
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
0 Comments
Accepted Answer
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.
4 Comments
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.
More Answers (0)
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!