Error using odearguments (line 95) SOL returns a vector of length 18, 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.
1 view (last 30 days)
Show older comments
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.
this is my code.
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.
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
hour=t1:1:t2;
costhetaz;
costheta_w;
costheta_g;
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
0 Comments
Answers (1)
Alan Stevens
on 20 Dec 2020
Your calculations of rad_w and rad_g produce vectors of size 8x1 that lead to tdot being larger than 3x1. Should these two (i.e. rad_w and rad_g) be the sum of the values in the vectors?
9 Comments
Alan Stevens
on 20 Dec 2020
As I noted before, where you currently have
for hour=t1:1:t2
replace this with
if t>=t1 & t<=t2
then work out lot etc for this specific time using the value of t.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!