equation system with time and space derivation

hello everyone!

i have an equation system for a solar panel formed by 4 equations 3 of them are depending only on time variation and the last one is depending on time and space variation .could you help me to implement the last equation on matlab .does ode45 works to solve this equation system .

dTfluid/dz i: the variation of temperature in the coil

best regards

Answer 1

[Alan Stevens]

You could use the pdepe function (try doc pdepe). Alternatively, explicitly discretise z and use ode45 to solve the time dependence for each section.

[badia rt]

thanks for your answer,

i thinks that pdepe doesn't work in my situation.i am agree with you about using Z descritization but i didn't find how to implement it in matlab .could you help me in that sense.

[Alan Stevens]

If your z goes from say, z = 0 to z = L, and you want to have N equispaced values of temperature, then set

h = L/(N-1), so each section has a length h.

You can then replace dTflui/dz by (Tfluid(i) - Tfluid(i-1))/h

You also need to replace Tfluid elsewhere by Tfluid(i), and if Tglass, Tpvt and Twax vary with z, then they need to be replaced by Tglass(i) etc.

You will then have a system of ordinary differential equations to solve.

[badia rt]

thanks a lot for replying

i write my code and change the z equation like you said but i think that there is some thing wrong in implementation.

function dy=Prog5(t,y)

global a2 a3 a4 b2 b3 b4 c1 c2 c3 c4 d1 d2 TCOIL z n alphaverre Uv Ta sigma epsilonv Averre rho Cp rhopv Cppv deltav APVT deltapv

dy(1)= 0.08.*sin(t./6000);

dy(2)=0.105+7.8*10^-6.*y(1)-a2.*y(2)-(a3.*y(2).^4)+a4.*y(3);

dy(3)=((1-alphaverre).*y(1)/(rhopv*Cppv*deltapv*APVT))-b2.*y(2)+b3.*y(3)+b4.*y(4);

dy(4)=c1+c2.*y(3)-c3.*y(4)-c4.*y(5,i);

% THE EQUATION DEPENDING ON TIME AND SPACE DERIVATION

dy(5)=d1.*(y(4)-y(5,i))-d2.*(y(5,i)-y(5,i-1))/h;

dy=[dy(1);dy(2);dy(3);dy(4);dy(5)];

end

[Alan Stevens]

Here's a version that uses a simple Euler integration technique for the time integration.

The constants are quite arbitrary, so the output from the coding as it stands will be nonsense. You will have to replace the constants with your own.

Note that because dTwax depends on Tfluid, then it, too, is indirectly a function of z. Similarly for the others.

 % Data - arbitrary values - replace with real ones
a2 = 10^-4; a3 = 10^-4; a4 = 1; 
alphaverre = 0.01; rhopv = 1; Cppv = 1; deltapv = 1; APTV = 1;
b2 = 10^-4; b3 = 1; b4 = 10^-4;
c1 = 1; c2 = 1; c3 = 10^-3; c4 = 10^-3;
d1 = 1; d2 = 10^-2;

% Spatial geometry
L =  2 ;      % Insert length 
N =  20 ;      % insert number of temperature nodes
h = L/(N-1); % length of section (ie distance between temperature nodes

% Timespan
tf = 1;        % insert final time
n  = 100;       % insert number of timesteps for which you want output
dt = tf/n;     % timestep

% Set aside storage
Tglass = zeros(N,n); 
Tpvt   = zeros(N,n);
Twax   = zeros(N,n);
Tfluid = zeros(N,n);
t = zeros(1,n);
z = zeros(1,N);

% Initialise temperatures
Tglass(1,1)  = 100;
Tpvt(1,1)    = 80;
Twax(1,1)    = 60;
Tfluid(1,1) = 20;

for i = 2:N   % space loop
    
   Tglass(i,1) = 100;
   Tpvt(i,1)   = 80;
   Twax(i,1)   = 60;
   Tfluid(i,1) = 20;
   z(i) = (i-1)*h;
   
   for j = 2:n  % time loop
       
       % Use a simple Euler integration for simplicity
       t(j) = (j-1)*dt;
       sinterm = 0.08*sin(t(j)/6000);
       
       dTglassdt = 0.105+7.8*10^-6*sinterm -a2*Tglass(i,j-1) ...
                   -a3*Tglass(i,j-1)^4 + a4*Tpvt(i,j-1);
       dTpvtdt = (1-alphaverre)*sinterm/(rhopv*Cppv*deltapv*APTV) ...
                  - b2*Tglass(i,j-1) + b3*Tpvt(i,j-1) + b4*Twax(i,j-1);
       dTwaxdt =  c1 + c2*Tpvt(i,j-1) - c3*Twax(i,j-1) - c4*Tfluid(i,j-1);
       dTfluiddz = (Tfluid(i,j-1)- Tfluid(i-1,j-1))/h;
       dTfluiddt = d1*(Twax(i,j-1) - Tfluid(i,j-1)) - d2*dTfluiddz;
       
       Tglass(i,j) = Tglass(i,j-1) + dt*dTglassdt;
       Tpvt(i,j)   = Tpvt(i,j-1) + dt*dTpvtdt;
       Twax(i,j)   = Twax(i,j-1) + dt*dTwaxdt;
       Tfluid(i,j) = Tfluid(i,j-1) + dt*dTfluiddt;
    end
    
end

subplot(2,1,1)
plot(t,Tfluid),grid
xlabel('time'),ylabel('fluid')
subplot(2,1,2)
plot(z,Tfluid),grid
xlabel('space'),ylabel('fluid')

You might want to modify this to use a better integration routine (Heun, RK4 etc), but, with luck, it gives you a starter for ten!

CHECK THE CODING CAREFULLY! I did it rather hurriedly and it wouldn't surprise me if it contained errors.

[badia rt]

thank you very much Alan Stevens.i replaced the data with the real value but the results didn't give me a suitable curves ,i run the simulation for 10 hours (means 36000 seconds because the time in sinterm is second).here is the code

% Data - arbitrary values - replace with real ones
alphaverre=0.03;Uv=0.15; Averre=0.8571;Ta=25;sigma=5.670367*10^-8;epsilonv=0.88;rho=3000;
Cp= 500;kverre=1;rhopv=2330;Cppv=702.9;deltav=0.003;kPVT=149;APVT =0.8571;deltapv=0.0002;
m=0.12;Cpf=4180;
Tfin=25;rhocire=930;Cpcire=2100;h=30;mf=0.12;cpf=4180;rhof=1000;Acoil=6.3*10^-5;d=2*sqrt(Acoil/pi);
hf=166.216;AWAX=0.8571;ewax=0.002;
a2=(((5.7+3.8.*Uv).*Averre)/ rho*Cp*Averre*deltav)+(kverre.* Averre)/( rho*Cp*Averre*deltav*deltav);
a3=(sigma*epsilonv* Averre)/( rho*Cp*Averre*deltav);
a4=(kverre* Averre)/(rho*Cp*Averre*deltav*deltav);
b2=(Averre*kverre)/((rhopv*Cppv*deltapv*APVT)*deltav);
b4=kPVT*APVT/(rhopv*Cppv*deltapv*APVT*deltapv);
c1=((h*APVT*Ta)+(m*Cpf* Tfin))/(rhocire*Cpcire*AWAX*ewax);
c2=(kPVT*APVT)/(rhocire*Cpcire*AWAX *ewax*deltapv);
c3=(kPVT*APVT)/((rhocire*Cpcire*AWAX *ewax)*deltapv)+h*APVT/(rhocire*Cpcire*AWAX *ewax);
c4=(m*Cpf)/(rhocire*Cpcire*AWAX *ewax);
d1=pi*d*hf/rhof*cpf*Acoil;
d2=(mf*cpf)/(cpf*rhof*Acoil);
b3= (kverre*Averre/((rhopv*Cppv*APVT*deltapv)*deltav))-(kPVT*APVT/((rhopv*Cppv*deltapv*APVT)*deltapv));
%The constant's value
%a2 = 0.075; a3 = 1.1*10^-11; a4 =0.075; 
%alphaverre = 0.03; rhopv = 2330; Cppv = 702.9; deltapv = 0.0002; APTV = 0.8571;
%b2 = 1.018; b3 = -2273.47; b4 = 2274.5;
%c1 = 3.78; c2 = 190.72; c3 = 190.72; c4 = 0.15;
%d1 = 0.00123; d2 = 1.9;
% Spatial geometry
L =  2 ;      % Insert length 
N =  20 ;      % insert number of temperature nodes
h = L/(N-1); % length of section (ie distance between temperature nodes
% Timespan
tf = 36000;        % insert final time
n  = 10;       % insert number of timesteps for which you want output
dt = tf/n;     % timestep
% Set aside storage
Tglass = zeros(N,n); 
Tpvt   = zeros(N,n);
Twax   = zeros(N,n);
Tfluid = zeros(N,n);
t = zeros(1,n);
z = zeros(1,N);
% Initialise temperatures
Tglass(1,1)  = 29;
Tpvt(1,1)    = 25;
Twax(1,1)    = 20;
Tfluid(1,1) = 19;
for i = 2:N   % space loop
    
   Tglass(i,1) = 29;
   Tpvt(i,1)   = 25;
   Twax(i,1)   = 20;
   Tfluid(i,1) = 19;
   z(i) = (i-1)*h;
   
   for j = 2:n  % time loop
       
       % Use a simple Euler integration for simplicity
       t(j) = (j-1)*dt;
       sinterm = 0.08*sin(t(j)/6000);
       
       dTglassdt = 0.105+7.8*10^-6*sinterm -a2*Tglass(i,j-1) ...
                   -a3*Tglass(i,j-1)^4 + a4*Tpvt(i,j-1);
       dTpvtdt = (1-alphaverre)*sinterm/(rhopv*Cppv*deltapv*APTV) ...
                  - b2*Tglass(i,j-1) + b3*Tpvt(i,j-1) + b4*Twax(i,j-1);
       dTwaxdt =  c1 + c2*Tpvt(i,j-1) - c3*Twax(i,j-1) - c4*Tfluid(i,j-1);
       dTfluiddz = (Tfluid(i,j-1)- Tfluid(i-1,j-1))/h;
       dTfluiddt = d1*(Twax(i,j-1) - Tfluid(i,j-1)) - d2*dTfluiddz;
       
       Tglass(i,j) = Tglass(i,j-1) + dt*dTglassdt;
       Tpvt(i,j)   = Tpvt(i,j-1) + dt*dTpvtdt;
       Twax(i,j)   = Twax(i,j-1) + dt*dTwaxdt;
       Tfluid(i,j) = Tfluid(i,j-1) + dt*dTfluiddt;
    end
    
end
figure(1)
plot (t,sinterm)
title('FLUX SOLAIRE');
xlabel('Time t');
ylabel('puissance solaire');
legend('G')
hold on
figure(2)
plot(t,Tglass,'r')
hold on
plot(t,Tpvt,'y')
hold on
plot(t,Twax,'g')
hold on
title('temperature variation in time');
xlabel('Time t');
ylabel('temperature T  °C');
legend('T_g','T_pvt','T_wax','T_fluide')
figure(3)
subplot(2,1,1)
plot(t,Tfluid),grid
xlabel('time'),ylabel('Tfluid')
subplot(2,1,2)
plot(z,Tfluid),grid
xlabel('space'),ylabel('Tfluid')

the results( sinterm,Tglass,Tpvt,Twax,Tfluid) must be like this curve

[Alan Stevens]

In order to get the above curve for G, the sinterm expression needs to be changed to

0.08*(sin(t(j)/tau-pi/2)+1)*950/0.16+50;

or something close to this. The result is then

However, this doesn't improve the remaining calculations, even with an increase in the number of sections and timesteps.

I've just realized that by assuming uniform temperatures along z for all times, there is nothing to make dTfluid/dz change! More information on what makes Tfluid vary with z is required.

[badia rt]

thanks a lot Alan Stevens what is the value of "tau" that make this curve for solar flux G .is it the same shape of the curve for the other output??means Tglass Tpvt Twax Tfluid

about dTfluid/dz i didn't have other equations but we can make an assumption that the temperature is uniforme in each section or there is a linear variation

[Alan Stevens]

Sorry, I should have noted that tau is just 6000, as in your original expression for sinterm.

There was no sensible solution for the other temperatures, but I didn't try further because of uncertainty in the z dependence. if the temperature is uniform along z then the dTfluid/dT will be zero (as it is in what I've supplied so far). If there is a inear variation then you will need to specity its magnitude.