using ode45 to solve 3 ODEs with solar radiation

3 次查看（过去 30 天）
Hasan Abbas
Hasan Abbas 2020-12-23
评论： Jan 2020-12-27
I am trying to solve 3 ODEs by using ode45. i dont know if im doing this right or not, I appreciate any small help.
my main problem is how to include the transient effect of solar radiation. I tried to solve it with constant solar radiation and it worked fine.
this is my code and when i run it i get this error.
Thanks in advance!
Error using odearguments (line 95)
SOL returns a vector of length 17, but the length of initial conditions vector is 3. The vector returned by SOL and
the initial conditions vector must have the same number of elements.
Error in ode45 (line 115)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in test (line 51)
t1=8; %start time
t2=15; %end time
tspan=[0,(t2-t1)*3600];
IC=[38,40,37]; %initial temperatures
[time,temp]=ode45(@sol,tspan,IC);
function tdot =sol(t,T)
t1=8;
t2=15;
vol=4;
mair=vol*1.225;
cp=1005;
Qh=0;%220;
Tac=15;
mac=0.09;
Uinwall=4;
Uoutwall=4.5;
Uinmass=12;
Amass=3;
Tamb=45;
cpwall=600;
cpmass=1450;
mwall=900;
mmass=50;
emiss=0.9;
Tsky=22.6;
abswall=0.19;
absmass=1;
trans=0.8;
A=[2.479,0 %roof
0.864,0.864 %windsheild
0.787,0.787 %rear window
1.533+0.224+0.266,0.224+0.266 %right side
1.533+0.224+0.266,0.224+0.266]; %left side
Aw=A(:,1); %area of walls
Ag=A(:,2); % area of glass
Awall=sum(Aw,'all');
AHwall=Aw(1,:);
lat=33.3152; %latitude
long=44.3661; % longitude
gmt=3;
mint=0;
tiltang_w=[0 % tilt angle of the walls
60
45
90
90];
zs_w=[0 % zenith angle of the walls
0
180
90
270];
tiltang_g=[60 % tilt angle of the glass
45
90
90];
zs_g=[0 % zenith angle of the glass
180
90
270];
costhetaz=[]; % cosine of the zenith angle
costheta_w=[]; %incident angle of solar radiation on the walls
costheta_g=[]; %incident angle of solar radiation on the glass
for hour=t1:1:t2
lot=hour+mint/60;%local time of the city
dn=190; % day of the year
a=1160+75*sind((dn-275)*360/365); %constants of the equation
b=0.174+0.035*sind((dn-100)*360/365); %constants of the equation
c=0.095+0.04*sind((dn-100*360/365)); %constants of the equation
delta=23.45*sind((dn+284)*360/365); % Declination angle
bn=(dn-1)*360/365;
eot=2.2918*(0.0075+0.1686*cosd(bn)-...
3.2077*sind(bn)-1.4615*cosd(2*bn)-...
4.089*sind(2*bn)); % equation of time
lst=15*gmt; % local standard time
st=(lot+(eot/60)+(long-lst)/15); % solar time
omega=15*(12-st); % hour angle
costhetaz_i=cosd(delta).*cosd(lat).*cosd(omega)+...
sind(delta).*sind(lat); % cosine of the zenith angle
costheta_wi=sind(delta).*sind(lat).*cosd(tiltang_w)-...
sind(delta).*cosd(lat).*sind(tiltang_w).*cosd(zs_w)+...
cosd(delta).*cosd(lat).*cosd(tiltang_w).*cosd(omega)+...
cosd(delta).*sind(lat).*sind(tiltang_w).*cosd(zs_w).*cosd(omega)+...
cosd(delta).*sind(tiltang_w).*sind(zs_w).*sind(omega); %incident angle of solar radiation on the walls
costheta_gi=sind(delta).*sind(lat).*cosd(tiltang_g)-...
sind(delta).*cosd(lat).*sind(tiltang_g).*cosd(zs_g)+...
cosd(delta).*cosd(lat).*cosd(tiltang_g).*cosd(omega)+...
cosd(delta).*sind(lat).*sind(tiltang_g).*cosd(zs_g).*cosd(omega)+...
cosd(delta).*sind(tiltang_g).*sind(zs_g).*sind(omega);%incident angle of solar radiation on the glass
costhetaz=[costhetaz;costhetaz_i];
costheta_w=[costheta_w;costheta_wi];
costheta_g=[costheta_g;costheta_gi];
end
minLength = min(length(costhetaz), length(costheta_w)); % i wrote this line becasue i didnt get the same number of values for costheatz,costheta_w and costheta_g
costhetaz = costhetaz(1:minLength);
costheta_w = costheta_w(1:minLength);
costheta_g = costheta_g(1:minLength);
[idir]=a.*exp(-b./costhetaz); % direct solar radiation
costheta_w(costheta_w<0)=0; % this line to neglect the negative values of costheta_w and costheta_g
costheta_g(costheta_g<0)=0;
rad_w=idir.*costheta_w; % solar radiation on the walls
rad_g=idir.*costheta_g; % solar radiation on the glass
totalrad_w=0;
totalrad_g=0;
for i=1:1:5
totalrad_w=totalrad_w+rad_w.*Aw(i,1);
end
for i=1:1:5
totalrad_g=totalrad_g+rad_g.*Ag(i,1);
end
Q_ac=0;%(mac*cp*(Tac-T(1)));
Qin_w=(Uinwall*Awall*(T(1)-T(2)));
Qin_m=(Uinmass*Amass*(T(1)-T(3)));
Qout_w=(Uoutwall*Awall*(Tamb-T(2)));
Qin_rad=(5.67*10^-8*Amass*((T(3)+273)^4-(T(2)+273)^4));
Qout_rad=-(5.67*10^-8*emiss*AHwall*((T(2)+273)^4-(Tsky+273)^4));
Qabssolar=abswall*totalrad_w;
Qtranssolar=trans*absmass*totalrad_g;
tdot=[(Qh+Q_ac-Qin_w-Qin_m)/(mair*cp)
(Qin_w+Qout_w+Qin_rad+Qabssolar-Qout_rad)/(mwall*cpwall)
(Qin_m+Qtranssolar-Qin_rad)/(mmass*cpmass)];
end
  5 个评论
显示 3更早的评论
Hasan Abbas
Hasan Abbas 2020-12-27
编辑：Hasan Abbas 2020-12-27
hello Jan
it is the same problem mentioned above, i solved the pervious error, i found that the terms "Qabssolar" and "Qtranssolar" were vectors of length (1*n) where n is the same as time span, so i searched for a solution and wrote this code and it worked well
Qabs_t = linspace(t1*3600,t2*3600,h);
Qtrans_t = linspace(t1*3600,t2*3600,h);
Qabssolar=abswall*totalrad_w;
Qtranssolar=trans*absmass*totalrad_g;
Qabssolar = interp1(Qabs_t, Qabssolar, t);
Qtranssolar = interp1(Qtrans_t, Qtranssolar, t);
i guess it`s right, because the result were as i expected. do you think it`s right or wrong?
this is my full code
t1=8;
t2=15;
h=((t2-t1)*3600)/60;
tspan=linspace(t1*3600,t2*3600,h);
IC=[32,33,35];
[time,temp]=ode45(@sol,tspan,IC);
function tdot =sol(t,T)
t1=8;
t2=15;
h=((t2-t1)*3600)/60;
n=((t2-t1)*3600)/h;
vol=4;
mair=vol*1.225;
cp=1005;
Qh=0;%220;
Tac=15;
mac=0.09;
Uinwall=4;
Uoutwall=4.5;
Uinmass=12;
Amass=3;
Tamb=46;
cpwall=600;
cpmass=1450;
mwall=900;
mmass=50;
emiss=0.9;
Tsky=22.6;
abswall=0.19;
absmass=1;
trans=0.8;
Aw=[2.479 %roof
0.864 %windsheild
0.787 %rear window
1.533+0.224+0.266 %right side
1.533+0.224+0.266]; %left side
Ag=[0.864 %windsheild
0.787 %rear window
0.224+0.266 %right windows
0.224+0.266];%left windows
Awall=sum(Aw,'all');
AHwall=Aw(1,:);
lat=33.3152;
long=44.3661;
gmt=3;
dn=212; % day of the year
a=1160+75*sind((dn-275)*360/365); %constants of the equation
b=0.174+0.035*sind((dn-100)*360/365); %constants of the equation
c=0.095+0.04*sind((dn-100*360/365)); %constants of the equation
delta=23.45*sind((dn+284)*360/365); % Declination angle
bn=(dn-1)*360/365;
mint=0;
tiltang_w=[0
60
45
90
90];
zs_w=[0
0
180
90
270];
tiltang_g=[60
45
90
90];
zs_g=[0
180
90
270];
costhetaz=[];
costheta_w=[];
costheta_g=[];
for hour=t1:1/n:t2-1/60
lot=hour+mint/60;%local time of the city
eot=2.2918*(0.0075+0.1686*cosd(bn)-...
3.2077*sind(bn)-1.4615*cosd(2*bn)-...
4.089*sind(2*bn)); % equation of time
lst=15*gmt; % local standard time
st=(lot+(eot/60)+(long-lst)/15); % solar time
omega=15*(12-st); % hour angle
costhetaz_i=cosd(delta).*cosd(lat).*cosd(omega)+...
sind(delta).*sind(lat); % cosine of the zenith angle
costheta_wi=sind(delta).*sind(lat).*cosd(tiltang_w)-...
sind(delta).*cosd(lat).*sind(tiltang_w).*cosd(zs_w)+...
cosd(delta).*cosd(lat).*cosd(tiltang_w).*cosd(omega)+...
cosd(delta).*sind(lat).*sind(tiltang_w).*cosd(zs_w).*cosd(omega)+...
cosd(delta).*sind(tiltang_w).*sind(zs_w).*sind(omega);
costheta_gi=sind(delta).*sind(lat).*cosd(tiltang_g)-...
sind(delta).*cosd(lat).*sind(tiltang_g).*cosd(zs_g)+...
cosd(delta).*cosd(lat).*cosd(tiltang_g).*cosd(omega)+...
cosd(delta).*sind(lat).*sind(tiltang_g).*cosd(zs_g).*cosd(omega)+...
cosd(delta).*sind(tiltang_g).*sind(zs_g).*sind(omega);
costhetaz=[costhetaz;costhetaz_i];
costheta_w=[costheta_w;costheta_wi];
costheta_g=[costheta_g;costheta_gi];
end
idir=a.*exp(-b./costhetaz); % direct solar radiation
idif=c.*idir;
i_tot=idir+idif;
costheta_w(costheta_w<0)=0;
costheta_g(costheta_g<0)=0;
rad_w=zeros(5,h);
for i=1:h
for j=1:5
rad_w(j,i)=rad_w(j,i)+ i_tot(i,1).*costheta_w(j,1);
end
end
rad_g=zeros(4,h);
for i=1:h
for j=1:4
rad_g(j,i)=rad_g(j,i)+ i_tot(i,1).*costheta_g(j,1);
end
end
totalrad_w=zeros(5,h);
totalrad_g=zeros(5,h);
for i=1:h
for j=1:1:5
totalrad_w(j,i)=totalrad_w(j,i)+rad_w(j,i).*Aw(j,1);
end
end
totalrad_w=sum(totalrad_w,1);
for i=1:h
for j=1:1:4
totalrad_g(j,i)=totalrad_g(j,i)+rad_g(j,i).*Ag(j,1);
end
end
totalrad_g=sum(totalrad_g,1);
Qabs_t = linspace(t1*3600,t2*3600,h);
Qtrans_t = linspace(t1*3600,t2*3600,h);
Qabssolar=abswall*totalrad_w;
Qtranssolar=trans*absmass*totalrad_g;
Qabssolar = interp1(Qabs_t, Qabssolar, t);
Qtranssolar = interp1(Qtrans_t, Qtranssolar, t);
Q_ac=0;%(mac*cp*(Tac-T(1)));
Qin_w=(Uinwall*Awall*(T(1)-T(2)));
Qin_m=(Uinmass*Amass*(T(1)-T(3)));
Qout_w=(Uoutwall*Awall*(Tamb-T(2)));
Qin_rad=(5.67e-8*Amass*((T(3)+273)^4-(T(2)+273)^4));
Qout_rad=(5.67e-8*emiss*AHwall*((T(2)+273)^4-(Tsky+273)^4));
tdot=[(Qh+Q_ac-Qin_w-Qin_m)/(mair*cp)
(Qin_w+Qout_w+Qin_rad+Qabssolar-Qout_rad)/(mwall*cpwall)
(Qin_m+Qtranssolar-Qin_rad)/(mmass*cpmass)];
end
Jan
Jan 2020-12-27
I do not have information about the purpose of your code. Therefore I cannot estimate, if the code is working as wanted or not. At least the variables Tac and mac are not used anywhere.
The expression "(T(3)+273)^4" seems to be very sensitive. What about using more digits like 273.15?
Qabssolar = interp1(Qabs_t, Qabssolar, t);
Qtranssolar = interp1(Qtrans_t, Qtranssolar, t);
Linear interpolations are not smooth. The stepsize control of ODE45 requires a smooth function. This can cause serious troubles, e.g. results, which are dominated by rounding errors. A stopping integration is possble also. Running a numerical integrator with not matching inputs is not a reliable method from a scientific point of view. If this is a part of a PhD thesis, I would reject it, if I am your professor.
The code does not allow to check the stability of the solution by varying the initial conditions and parameters. Did you validate the applied model?

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Mischa Kim
Mischa Kim 2020-12-24
Hi Hasan,
based on the error message I suggest looking into tdot at the bottom of your ode function. Apparently, this is a 17x1 vector and by browsing through your code it seems that Qabssolar and Qtranssolar are also vectors, not scalars.
I suggest you set a breakpoint next to the last end command in your code above and run the code. The execution of the code will halt when the breakpoint is reached and you can hover over the variables to see value and size. Once you have confirmed that Qabssolar and Qtranssolar are indeed the cause of the issue you can work your way back up and isolate the source of the problem.
  1 个评论
Hasan Abbas
Hasan Abbas 2020-12-24
编辑：Hasan Abbas 2020-12-24
Thanks for your answer, i checked it and it is indeed a vector. but i have another problem, i will appreciateit if you can help me. I want the values of rad_w and rad_g to change by taking values from t1 to t2 and i want the ode45 to evaluate the 'sol' function at each time span corresponding to that time. so i have time span from t1 to t2 and i have for loop from t1 to t2.
Also, in that for loop i get vector of 8*1 for costhetaz but i get a vector of 40*1 for costheta_w and a vector of 32*1 for costheta_g because i have vectors of 5*1 for the walls and 4*1 for the glass, how can i make it all simultaneous so at each time span it take 1 value of costheaz and 5 values of costheta_w and 4 values of costheat_g?
if i have done this wrong, please explain it to me
thanks again

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by