Error using pdepe: Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative
28 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Lorenz.zz
le 10 Juil 2024 à 15:00
Modifié(e) : Torsten
le 12 Juil 2024 à 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:
![](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 commentaires
Réponse acceptée
Torsten
le 10 Juil 2024 à 15:23
Modifié(e) : Torsten
le 10 Juil 2024 à 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 commentaires
Torsten
le 12 Juil 2024 à 9:30
Modifié(e) : Torsten
le 12 Juil 2024 à 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.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Eigenvalue Problems dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!