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

14 次查看（过去 30 天）
Widitha Samrakoon
Widitha Samrakoon 2016-2-18
评论： Torsten 2016-2-22
I get this error in writing the code.... The code I wrote is as follows
function pdeop global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a; m = 0; x = linspace(0,1.5,31); t = linspace(0,3600,3601); sol = pdepe(m,@pdeoppde,@pdeopic,@pdeopbc,x,t); u1 = sol(:,:,1); u2 = sol(:,:,2); u3=sol(:,:,3); u4=sol(:,:,4);
% -------------------------------------------------------------- function [c,f,s] = pdeoppde(x,t,u,DuDx) c = [A*void*(1-hl)*Rhog;A*void*hl*Rhol;A*void*(1-hl)*Rhog;A*void*hl*Rhol]; f = [-G;L;-G;L] .*u; Cs=Rhog*Kos*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff); Ct=Rhog*Kot*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff); ysstar=-4.2913*u(2)^4+11.791*u(2)^3-12.01*u(2)^2+5.508*u(2); ytstar=0.0023*exp(6.0058*u(4)); M=Cs*(ysstar-u(1)); N=Ct*(ytstar-u(3)); s = [M;-M;N;-N]; end % -------------------------------------------------------------- function u0 = pdeopic(x) u0 = [0;0.5;0;0.5]; end % -------------------------------------------------------------- function[pl,ql,pr,qr] = pdeopbc(xl,ul,xr,ur,t) pl = [ul(1);0;ul(3);0]; ql = [0;0;0;0]; pr = [ur(1)-0.5;0;ur(3)-0.5;0]; qr = [0;0;0;0]; end end
Can anyone please suggest me a solution?
  3 个评论
显示 1更早的评论
Torsten
Torsten 2016-2-19
... and best include the PDE system with boundary and initial conditions in mathematical notation as pdf document such that we are able to compare it to your code.
Best wishes
Torsten.
Widitha Samrakoon
Widitha Samrakoon 2016-2-19
编辑：Widitha Samrakoon 2016-2-19
My code contains two functions named paracal.m which calculated parameters to be used in the pde system and then the main function named pdeop.m. First paracal.m is run and then pdeop.m is run. First I will paste the code of paracal.m and then pdeop.m. The error occurs when running pdeop.m.
function paracalc
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a;
A=pi*0.0254^2*0.25;
a=1710;
L=7.2888*10^-5;
G=2.4296*10^-3;
Rhog=842.77;
Rhol=875;
void=0.86;
Usg=G/(Rhog*A);
Usl=L/(Rhol*A);
hl=(1.5688*Usl^(1/3))*(1+4.1435*Usg^0.05)*(1-exp(-88.9560*Usl^(0.4)))^(2/3);
Kos=(1278102.274+6643.8785*Usg^-0.5)^-1;
Kot=(1620456.645+10029.7854*Usg^-0.5)^-1;
aeff=a*(1-exp(-88.9560*Usl^0.4));
Ugeff=1.6444*Usg/(1-hl);
Uleff=1.6444*Usl/hl;
end
Now I will paste the code of pdeop.m which gives the error
function pdeop
global A void G L Usg Usl Rhog Rhol Kos Kot aeff Ugeff Uleff hl a;
m = 0;
x = linspace(0,1.5,31);
t = linspace(0,3600,61);
sol = pdepe(m,@pdeoppde,@pdeopic,@pdeopbc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);
u3=sol(:,:,3);
u4=sol(:,:,4);
function [c,f,s] = pdeoppde(x,t,u,DuDx)
c = [A*void*(1-hl)*Rhog;A*void*hl*Rhol;A*void*(1-hl)*Rhog;A*void*hl*Rhol];
f =[-G;L;-G;L].*u+[-0.001;0.001;-0.001;0.001].*DuDx;
Cs=Rhog*Kos*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff);
Ct=Rhog*Kot*aeff*A*sqrt(Rhog/(Rhol-Rhog))*(Ugeff/Uleff);
ysstar=-4.2913*u(2)^4+11.791*u(2)^3-12.01*u(2)^2+5.508*u(2);
ytstar=0.0023*exp(6.0058*u(4));
M=Cs*(ysstar-u(1));
N=Ct*(ytstar-u(3));
s = [M;-M;N;-N];
function u0 = pdeopic(x)
u0 = [0;0.5;0;0.5];
end
function[pl,ql,pr,qr] = pdeopbc(xl,ul,xr,ur,t)
pl = [ul(1);0;ul(3);0];
ql = [0;0;0;0];
pr = [0;ur(2)-0.5;0;ur(4)-0.5];
qr = [0;0;0;0];
end
end
Please note that my system of PDEs has four elliptical equations but the problem is all four of them do not have second order spatial derivatives. Therefore in my system, the flux vector doesn't have a DuDx term but U term only. By a search in the internet I got to know that pdepe function cannot work with such equations.However in order to get rid of this, as suggested in an answer to a similar question, I included a second order spatial derivative to all four equations which you can see as a separate column vector in the flux vector as [-0.001;0.001,-0.001,0.001]. I hoped this would solve the problem. But the error message keeps popping out. I am really helpless. Please help.
I am attaching a pdf file of my PDE system with the initial conditions and boundary conditions.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Bill Greene
Bill Greene 2016-2-20
I see two problems.
First, your boundary conditions are incorrectly defined. They should be:
pl = [ul(1);ul(2);ul(3);ul(4)];
pr = [ur(1);ur(2)-0.5;ur(3);ur(4)-0.5];
(ql and qr are correct). All four of your equations are hyperbolic which pdepe is not designed to handle. But, you have added some artificial diffusion (the small second derivative terms). That is a good approach but because of the negative signs on some of the coefficients you have de-stabilized the system instead stabilizing it. Try this for the vector of diffusion coefficients, instead:
eps=.001;
[eps eps eps eps]'
  2 个评论
Widitha Samrakoon
Widitha Samrakoon 2016-2-21
Thank you Bill for the hint. However, this system only has four know boundary conditions ( which should be enough to solve since the order of spatial derivatives is one in each). ul(2), ul(4) and ur(1) and ur(3) are not known. that's why I put 0 for both p and q in place of those variable. The objective of solving this system is finding the variation of those values with time....However, I will adjust my diffusion derivative and check and reach you back. Thank you
Bill Greene
Bill Greene 2016-2-21
OK, I took a guess at what your intention was with respect to the boundary conditions. pdepe requires boundary conditions at both ends and p=q=0 is not acceptable. I suggest p=0 and q=1 at the ends where you really don't want a BC. This will set your "artificial" flux to zero there but have minimal impact on your "real" PDE.

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2016-2-22
As Bill already pointed out, pdepe is not designed to solve PDEs without 2nd derivative terms.
You will have to discretize the spatial first derivatives (e.g. by a first-order upwind scheme) and solve the resulting system of ordinary differential equations using ODE15S, e.g. . Look up "method of lines" for more details.
Best wishes
Torsten.

类别

在 Help Center 和 File Exchange 中查找有关 Geometry and Mesh 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by