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];
endI 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 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 pde, pdepe MATLAB Answers — New Questions
