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):
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:
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.
Any idea or request for further explanations will be welcome :) :) :)
Thank you in advance
Giulio
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?
