here is a Matlab beginner banging his head on the wall...

I am trying to solve a system of partial differential equations in Matlab, with both derivatives in time and space domains. I am using the pdepe function for that.

The system is, to be simple, a sort of solar thermal panel, made of three layers: an absorber plate, a fluid layer of running water and a back insulation layer. each of these layers is represented by a differential equation. Of course, many assumptions are implicit in the model, e.g. the model is here considered 1D (no edge effects...), the temperature disuniformity inside each layer is neglected ecc.

The differential equations that describe the system are (at the present level of complexity): enter image description here

where CXX are the thermal capacitances and hXX are the heat transfer coefficients between the different layers of the thermal model. Vfl is the fluid velocity, Tia is a known ambient temperature and Psolar is a heating power due to solar irradiance.

As initial conditions: the temperature of the recirculating water varies linearly from an inlet temperature to a first-attempt outlet temperature of 12°C

the temperature of the solar absorber varies linearly from 16.5 to 19 (respectively, averages between the inlet temperature of the fluid and the ambient temperature and between the outlet temperature of the fluid and the ambient temperature)

analogously, the temperature of the back insulation layer varies linearly from 16.5 to 19.

The spatial boundary conditions are: enter image description here

where Tfl,inlet and Tfl,outlet are the inlet and outlet temperatures of the circulating water.

In Matlab this becomes (not all code lines are reported)

function SolarAbsorber_TimeSpace

m = 0;

x = linspace(0,2.5,30);

t = linspace(0,7200,100); 

sol = pdepe(m,@pd_SolAbs_pde,@pde_SolAbs_ic,@pde_SolAbs_bc,x,t);

T_met = sol(:,:,1); 

T_w = sol(:,:,2);

T_is = sol(:,:,3);

function [c,f,s] = pd_SolAbs_pde(x,t,u,DuDx)
% I do  not report here the initialization of the coefficients

c = [C_met;C_w;C_is];                         

f = [ 0; 0; 0];

s1 = (h_ai_met)*(T_ai - u(1))+ h_met_w*(u(2) - u(1))+ Gv*tau*abs);

s2 = (h_met_w*(u(1) - u(2))+ h_is_w*(u(3) - u(2)))- DuDx(2)*vx;

s3 = (h_is_w*(u(2) - u(3)) + h_is_ai*(T_ai - u(3)));

s = [s1; s2; s3];                                

% --------------------------------------------------------------------------  
function u0 = pde_SolAbs_ic(x)

alfa =((12-7)/2.5);

beta = ((26+12)/2 -(26+7)/2)/2.5;

T_met_0 = beta*x+(26+7)/2; %°C

T_w_0 = alfa*x+7;    %°C

T_is_0 = beta*x+(26+7)/2; %°C

u0 = [T_met_0; T_w_0; T_is_0];

% --------------------------------------------------------------------------  
function [pl,ql,pr,qr] = pde_SolAbs_bc(xl,ul,xr,ur,t,u0,u)

pl = [ul(1)- (26+7)/2; ul(2)-7; ul(3)-(26+7)/2];

ql = [0;0;0];  

pr = [ur(1)-(26+12)/2; ur(2)-12; ur(3)-(26+12)/2];

qr = [0;0;0];

I'd like to describe the profile of the water temperature, giving a constant mass flow rate, a constant inlet temperature, variable forcing boundary conditions such as solar radiation, the initial conditions for the three layers etc.

In the first simulation, I just assume that solar radiation is constant and, starting at time t=0, it heats up the absorber while the circulating water removes the heat.

Here my questions:

1 Matlab_pdepe asks me to define boundary conditions also on the "right side" (outlet of the solar panel), which is actually unknown!

2 The graph shows that in correspondence of the boundaries, the profile gets unstable. This does not change significantly with the mesh size nor with the integration time span. Please consider that the left boundary condition of the absorber plate is fixed by the heat balance with the water inlet temperature and the ambient temperature, both supposed constant.

enter image description here

Any idea or request for further explanations will be welcome :) :) :)

Thank you in advance


P.S. Thank you Bill Greene for the comment about the Peclet number and the issue about diffusion vs convection mechanisms. Actually, the water velocity is fixed (dependent on the fixed water mass flow). I have refined the mesh in order to decrease the Peclet number and the result is that the profile tends to stability, but for the regions in proximity of the boundaries. Of course, now the point is understanding more in depth what happened ;)

However, the results are not physically consistent: the water passes from 7 to 11.7°C in less than 1cm. In a similar way, close to the outlet, the temperature oscillates to fit with the ooutlet boundary condition :(

I still would like to fix only the inlet temperature and the water velocity, letting the outlet temperature free to vary according to the physics of the system. Any suggestions about the "right" boundary conditions for the three layers?

enter image description here enter image description here

  • 2
    $\begingroup$ Welcome to Scicomp.SE! It's difficult to know what's going on without actually seeing (in mathematical notation) the PDE you are trying to solve (including all boundary and initial conditions). If (1) is your main point: you need a boundary condition at every side, otherwise the problem does not make sense mathematically. However, there are different conditions modeling different effects (such as fixed flux instead of fixed temperature), so you should also give a bit more details about the model problem. $\endgroup$ Commented Apr 21, 2015 at 13:43
  • $\begingroup$ Please put the equations in the question, using latex -- not in a comment. $\endgroup$ Commented Apr 22, 2015 at 11:29
  • $\begingroup$ Unfortunately, I do not know the answer to your questions. However, I find you PDE model interesting. Did you create it, or do you perhaps have a literature/web reference? Thanks. $\endgroup$
    – user21
    Commented Apr 24, 2015 at 11:57
  • $\begingroup$ I based the model on some previous literature studies. See for example: Kamminga 1985 Int.Heat Mass Transfer Schnieders 1997 Solar Energy My interest is not on a very detailed description of temperature distribution (e.g. edge effects). I would use the model to predict the output temperture... but Matlab wants me to set the output temperature before simulating it :( Probably there's an easy way to solve this limit... $\endgroup$
    – Giulio
    Commented Apr 27, 2015 at 13:17

1 Answer 1


This looks like a classic case of the convective term in your equations being large relative to the diffusive (conductive) term. The pdepe function is basically designed to handle the opposite case-- where diffusion dominates. You can quantify the relative effects of convection and diffusion in your problem by computing the mesh Peclet number. This has been discussed in some detail in another question previously asked on this site: Does the finite element method impose any restrictions on the Peclet number for numerical stability?

  • $\begingroup$ Thanks to everyone who has answered so far :) This forum is really something!!! I have added some considerations and graphs in the question box. $\endgroup$
    – Giulio
    Commented Apr 22, 2015 at 13:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.