How to solve a system of parabolic partial differential equations in time and 1 D space

Hi,

I cannot seem to solve the PDF in the figures below. This is the code I have but it produces the wrong graphs:

function New
m=0;
x=0:0.5:15;
t=1:100:28800;

sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);

u1=sol(:,:,1);
u2=sol(:,:,2);

surf(x,t,u1)
title('u1(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')

function [c,f,s]=pdefun(x,t,u,DuDx)
P=1200;
R=0.208;
cv=0.313;
rho=7.5;
v=0.3;
A=70;
mdot=1.2;
k=0.016/1000;
alph=0.5;
rhoG=1600;
cG=0.7;
kG=0.7/1000;

c=[((P*cv/R)/u(1));(rhoG*cG)];

f=[((k*DuDx(1))-((P*cv/R)*(mdot*R/(P*A))));(kG*DuDx(2))];

s=[(u(2)-u(1));(u(1)-u(2))]*alph;

function u0=icfun(x)
u0=[300;300];

function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
pl = [ul(1)-773; 0]; 
ql = [0; 1]; 
pr = [0; 0]; 
qr = [1; 1];

I think my boundary conditions are incorrect but I don't know how to fit them.

I would appreciate any help with this problem, thank you in advance.

<<

>>

%

Answer 1

[Torsten]

One mistake in your settings:

f = [k*DuDx(1) ; kG*DuDx(2)];
s = [(u(2)-u(1))*alph-cv*mdot/A*DuDx(1);(u(1)-u(2))*alph];

And maybe in your boundary condition for T_A at x=0, you will need a ramp that linearly increases the temperature from 300 K to 773 K in, say, one minute.

Best wishes

Torsten.

[Lucy]

Hi Torsten,

Thank you very much for your answer I did not notice that mistake. When I make that change I still get funny looking graphs of under 300K and over 800K temperature when it should show a temperature change between 300 and 773. The graphs are also oscillating and I was reading that this could be because only a very few of the possible odeset() options are passed and to get a much more stable solution from the underlying ode15s solver the 'MaxOrder' option from the Default value of {5} to 2. I tried to this this but just got errors. Do you happen to know any other reason this might be happening?

Thank you again for your help.

Lucy

[Torsten]

What happens if you set "mdot" to a very small value ?

And did you try to use the ramp ? A discontinuous initial condition usually produces oscillations in the solution.

Best wishes

Torsten.

[Lucy]

Im sorry I am very new to Matlab, how would I add a ramp to the initial condition? All the examples I saw just have one line under the function:

function u0=icfun(x) u0=[300;300];

A smaller mdot creates more oscillations and negative temperatures.

[Lucy]

I got it to work! Thank you so much for finding that mistake for me!

function New
m=0;
x=0:.1:15;
t=0.1:200:28800;

options = odeset('MaxOrder',2);
sol = pdepe_opt(m,@pdefun,@icfun,@bcfun,x,t,options);

uA=sol(:,:,1);
uG=sol(:,:,2);

surf(x,t,uA)
title('uA(x,t)')
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,uG)
title('uG(x,t)')
xlabel('Distance x')
ylabel('Time t')

function [c,f,s]=pdefun(x,t,u,DuDx)
P=1200;
R=0.208;
cA=0.313;
A=70;
mdot=240;
kA=0.016/1000;
alph=0.5;
rhoG=1600;
cG=0.7;
kG=0.7/1000;

c=[((P*cA)/(R*u(1)));(rhoG*cG)];

f = [kA*DuDx(1);kG*DuDx(2)];

s = [(u(2)-u(1))*alph-cA*mdot/A*DuDx(1);(u(1)-u(2))*alph];

function u=icfun(x)
u=[300;300];

function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)

if t==0
        pl = [ul(1)-300; 0]; 
else 
        pl = [ul(1)-773; 0];
end

ql = [0; 1]; 
pr = [0; 0]; 
qr = [1; 1];

[Torsten]

Replace your function "bcfun" with

function[pl,ql,pr,qr]=bcfun(xl,ul,xr,ur,t)
a = 1000;
pl = [ul(1)-(300+(773-300)*min(t,a)/a) ; 0]; 
ql = [0; 1]; 
pr = [0; 0]; 
qr = [1; 1];

and leave everything else as it is.

If this does not help, "pdepe" cannot handle the discontinuity at x=0.

You will have to discretize your PDE system on your own and solve the resulting system of ordinary differential equations using ODE15S. Look up "method of lines" for more details.

Best wishes

Torsten.

[Torsten]

Increasing the number of discretization points together with the ramp for a=100 solved the problem:

x = linspace(0,15,3000);

Best wishes

Torsten.

[Lucy]

Hi Torsten,

Thank you so much for your help. It is now working how I believe it should.

Lucy