is it possible to solve 2 PDEs simultaneously using finite difference method in matlab

I want to model a thermal energy storage system where the temperatures of the HTF and the PCM are given as:

dTf/dt = Af*d^2Tf/dx^2 - Bf*dTf/dx - Cf(Tf - Tpcm) and

dTs/dt = Apcm*d^2Ts/dx^2 - Cpcm(Tpcm - Tf)

where Af, Apcm, Bf, Cf and Cpcm are constants.

Initial conditions are:

Tf(x,0) = 180 degC and Tpcm(x,0) = 25 degC

Boundary conditions:

At x=o: dTf/dx = 0 and Tpcm = 180 degC

At x=0.3: dTpcm/dx = 0 and Tf = 150 degC

I have used the pdepe, but I not happy with the temperature distributions i am getting, maybe I need to adjust my initial and boundary conditions. That's why I am trying to solve the 2 PDEs using finite difference method

[Torsten]

This is a typical problem for "pdepe". If you are not happy with its solution, you won't be happy with your finite-difference solution, either.

Best you post your "pdepe" code so that we can look over it.

[luckywell seyitini]

Thank you Torsten for the advice. Actually, I am still new to Matlab and i want to model a hybrid latent-sensible thermal storage system. The code I am sharing is for the latent section only and i intend to use similar PDEs for the sensible section of the storage system. The two shall be linked by the condition that, outlet temperature of HTF from the latent section is equal to inlet temperature of HTF into the sensible section.

I have used the following pdepe code, please check and share suggested adjustments:

pdex4()

Error using pdepe
Unexpected output of PDEFUN. For this problem PDEFUN must return three column vectors of length 2.

Error in solution>pdex4 (line 7)
sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);

function pdex4
m = 0;
x = [0:0.075:0.3];
t = [0:5:30];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

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

figure
plot(x,u1(end, :))
title ('Variation of HTF temperature')
xlabel('Distance x')
ylabel('Temperature T')

figure
plot(x,u2(end, :))
title('Variation of PCM temperature')
xlabel('Distance x')
ylabel('Temperature T')

figure 
imagesc(x,t,u1); colormap hot; colorbar; grid on
title('Temperature distribution in HTF');
xlabel('Distance x')
ylabel('Time t')

figure 
imagesc(x,t,u2); colormap hot; colorbar; grid on
title ('Temperature distribution in PCM');
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
end
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [1; 1]; 
Af = 0.000000007*u(1)^2-0.000000009*u(1)+0.0001;
Apcm = 0.000000000001*u(2)^2-0.0000000003*u(2)+0.00000008;
f = [Af; Apcm] .*DuDx ;
y = u(1) - u(2);
Cf = -0.000009*u(1)+0.0181;
Cpcm = 0.0000002*u(2)^2-0.00005*u(2)+0.0209;
Bf = -0.0000000002*u(1)^2+0.0000002*u(1)+0.0003;
F = Cpcm*y;
G=-Bf*DuDx-Cf*y;
s = [G; F] ;
end
% --------------------------------------------------------------
function u0 = pdex4ic(x);
u0 = [170;60];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [0; ul(2)-150]; 
ql = [1; 0]; 
pr = [ur(1)-150; 0]; 
qr = [0; 1];
end

[Torsten]

If you have an inflow at x = 0 for Tf (-Bf*DuDx(1)) , it's necessary that you also set the inflow temperature at this boundary. But you set dTf/dx = 0 at x = 0.

pdex4()

function pdex4
m = 0;
x = linspace(0,0.3,100);%[0:0.075:0.3];
t = [0:5:30];

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t);
u1 = sol(:,:,1);
u2 = sol(:,:,2);

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

figure
plot(x,u1(end, :))
title ('Variation of HTF temperature')
xlabel('Distance x')
ylabel('Temperature T')

figure
plot(x,u2(end, :))
title('Variation of PCM temperature')
xlabel('Distance x')
ylabel('Temperature T')

figure 
imagesc(x,t,u1); colormap hot; colorbar; grid on
title('Temperature distribution in HTF');
xlabel('Distance x')
ylabel('Time t')

figure 
imagesc(x,t,u2); colormap hot; colorbar; grid on
title ('Temperature distribution in PCM');
xlabel('Distance x')
ylabel('Time t')

figure
surf(x,t,u2)
title('u2(x,t)')
xlabel('Distance x')
ylabel('Time t')
end
% --------------------------------------------------------------
function [c,f,s] = pdex4pde(x,t,u,DuDx)
c = [1; 1]; 
Af = 0.000000007*u(1)^2-0.000000009*u(1)+0.0001;
Apcm = 0.000000000001*u(2)^2-0.0000000003*u(2)+0.00000008;
f = [Af; Apcm] .*DuDx ;
y = u(1) - u(2);
Cf = -0.000009*u(1)+0.0181;
Cpcm = 0.0000002*u(2)^2-0.00005*u(2)+0.0209;
Bf = -0.0000000002*u(1)^2+0.0000002*u(1)+0.0003;
F = Cpcm*y;
G=-Bf*DuDx(1)-Cf*y;
s = [G; F] ;
end
% --------------------------------------------------------------
function u0 = pdex4ic(x);
u0 = [170;60];
end
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t)
pl = [0; ul(2)-150]; 
ql = [1; 0]; 
pr = [ur(1)-150; 0]; 
qr = [0; 1];
end

[luckywell seyitini]

Thank you very much Torsten. Yes, there is an inflow of HTF at x = 0 and the assumption is the inlet temperature of the HTF is constant at 170 degC, that why I decided to set dTf/dx = 0 at x=0.

Another anormaly I am observing is, when time is increased from 30sec at some point, the Tf and Tpcm will go beyong the specified values in the initial and boundary conditions i.e it can get to negative values. What can be the problem?

[Torsten]

Yes, there is an inflow of HTF at x = 0 and the assumption is the inlet temperature of the HTF is constant at 170 degC, that why I decided to set dTf/dx = 0 at x=0.

If HTF is 170 degC at x=0 and flows into the domain, why don't you set

pl(1) = ul(1)-170; 
ql(1) = 0;

?

And are you sure you want to have

f = [Af; Apcm] .*DuDx ;

?

Remember that Af and Apcm depend on u. Thus d/dx(f*du/dx) does not equal f*d^2u/dx^2.

[luckywell seyitini]

Thank you for the suggestion on the boundary condition, I will try it.

Actually Af and Apcm are thermal diffusivity for the HTF and PCM respectively. The f, I have used to denote diffusivity for the fluid may be confusing. I will cross check, but I think it should be correct as it is. Thank you.

[Torsten]

I will cross check, but I think it should be correct as it is.

If the equations from your above question are correct, it's not correct as is.

You define

f = [Af; Apcm] .*DuDx ;

This means that the thermal diffusivity terms in MATLAB are

d/dx(Af*dTf/dx) = Af*d^2Tf/dx^2 + d/dx(Af) * dTf/dx
d/dx(Apcm*dTs/dx) = Apcm*d^2Ts/dx^2 + d/dx(Apcm) * dTs/dx

and the terms d/dx(Af) and d/dx(Apcm) are not zero because they depend on Tf and Ts.

But you claim your equations read

dTf/dt = Af*d^2Tf/dx^2 - Bf*dTf/dx - Cf(Tf - Ts) and 
dTs/dt = Apcm*d^2Ts/dx^2 - Cpcm(Ts - Tf),

not

dTf/dt = Af*d^2Tf/dx^2 + d/dx(Af) * dTf/dx - Bf*dTf/dx - Cf(Tf - Ts) and 
dTs/dt = Apcm*d^2Ts/dx^2 + d/dx(Apcm) * dTs/dx - Cpcm(Ts - Tf)

[luckywell seyitini]

okay, I get you Torsten, let me address that. thank you