下面的程序，用for循环计算t从9到16，每次到t=9.5左右就自己结束循环了。不知道为什么
clear;clc;
%初始温度
I=800;
T_sun=6000;
T_amb=20;%
T_g0=T_amb;
T_pv0=T_amb;
T_a=T_amb;
T_sky=T_amb-6;
T_g=T_g0;
T_pv=T_pv0;
T_a_in=T_amb;
T_a_out=2*T_a-T_a_in;
%结构尺寸
A_in=0.0238;
A_out=0.024;
delta_a=0.005;
D=0.5;
H=1.0;
A_r=delta_a*D;
A=D*H;
u_amb=1;%
delta_g=0.005;
epsilon_g=0.9;
alpha_g=0.1;
tau_g=0.9;
rho_g=2500;
c_g=840;
m_g=rho_g*delta_g*A;
lambda_a=0.026;
nu_a=1.58*10^(-5);
rho_a=1.18;
c_a=1100;
a=lambda_a/(rho_a*c_a);
Pr=nu_a/a;
m_a=rho_a*delta_a*A;
alpha_pv=0.9;
epsilon_pv=0.94;
eta_r=0.1;
rho_pv=2500;
c_pv=840;
delta_pv=0.005;
m_pv=rho_pv*delta_pv*A;
lambda_w=0.035;
delta_w=0.3;
n=0;
x=0;
delta_t=1/6480;%5秒
for t=9.5:delta_t:16
%环境参数拟合
T_a_in=T_amb;
u_amb=1;
I=359.9+19.13*cos(t*0.4052)-395.5*sin(t*0.4052);
%换热系数计算
h_amb_g=5.7+3.8*u_amb;
sigma=5.67*10^(-8);
g=9.80;
beta_g=1/((T_a+T_g)/2+273);
Ra_g=g*beta_g*(abs(T_a-T_g))*H^3/(nu_a*a);
if Ra_g<10^9
Nu_g=0.68+0.670*Ra_g^(1/4)/(1+(0.492/Pr)^(9/16))^(4/9);
else
Nu_g=(0.825+0.387*Ra_g^(1/4)/((1+(0.492/Pr)^(9/16))^(8/27)))^2;
end
h_a_g=Nu_g*lambda_a/H;
if isreal(h_a_g)==0
a=0;
end
h_pv_g=sigma/(1/epsilon_g+1/epsilon_pv-1);
beta_pv=1/((T_a+T_pv)/2+273);
Ra_pv=g*beta_pv*(abs(T_a-T_pv))*H^3/(nu_a*a);
%if Ra_pv<10^9
Nu_pv=0.68+0.670*Ra_pv^(1/4)/(1+(0.492/Pr)^(9/16))^(4/9);
%else
% Nu_pv=(0.825+0.387*Ra_pv^(1/4)/((1+(0.492/Pr)^(9/16))^(8/27)))^2;
% (0.825+0.387*Ra_pv^(1/4)/(1+(0.492/Pr)^(9/16))^(8/27))^2
%end
h_pv_a=Nu_pv*lambda_a/H;
h_g_a=h_a_g;
%自然对流
f_in=1.0;
f_out=1.5;
d=4*A_r/(2*delta_a+2*D);
Re=0.2*d/nu_a;
f=0.3164*Re^-0.25;
beta=1/(T_a+273);
u_a=sqrt(0.5*g*beta*(T_a_out-T_a_in)*H/(f_in*(A_r/A_in)^2+f_out*(A_r/A_out)^2+f*(H/d)));
u_a
h_g_pv=h_pv_g;
h_a_pv=h_pv_a;
h_w_amb=29;
R=delta_w/lambda_w+1/h_w_amb;
T_r=298.15;
q_amb_g=h_amb_g*A*(T_amb-T_g)
q_sky_g=sigma*epsilon_g*A*((T_sky+273)^4-(T_g+273)^4)
q_a_g=h_a_g*A*(T_a-T_g)
q_pv_g=h_pv_g*A*((T_pv+273)^4-(T_g+273)^4)
q_g_a=h_g_a*A*(T_g-T_a)
q_pv_a=h_pv_a*A*(T_pv-T_a)
q_a=rho_a*u_a*c_a*A_r*(T_a_out-T_a_in)
q_g_pv=h_g_pv*A*((T_g+273)^4-(T_pv+273)^4)
q_a_pv=h_a_pv*A*(T_a-T_pv)
q_amb_pv=(T_amb-T_pv)*A/R
E_pv=tau_g*A*I*eta_r*(1-0.0045*(T_pv+273.15-T_r))
%温度计算
T_1=(q_amb_g+q_sky_g+q_a_g+q_pv_g+alpha_g*A*I)*(delta_t*3600)/(m_g*c_g)
T_2=(q_g_a+q_pv_a-q_a)*(delta_t*3600)/(m_a*c_a)
T_3=(q_g_pv+q_a_pv+q_amb_pv+tau_g*alpha_pv*A*I-E_pv)*(delta_t*3600)/(m_pv*c_pv)
if T_2>1
break
end
T_g=T_g+T_1;
T_g=vpa(T_g,7);
T_a=T_a+T_2;
T_a=vpa(T_a,7);
T_pv=T_pv+T_3;
T_pv=vpa(T_pv,7);
T_a_out=2*T_a-T_a_in
n=n+1;
if mod(n-1,6480)==0;
x=x+1;
T_a_dsj(x)=[T_a];
T_g_dsj(x)=[T_g];
T_pv_dsj(x)=[T_pv];
T_a_out_dsj(x)=[T_a_out];
I_dsj(x)=[I];
data(1,x)=I_dsj(x);
data(2,x)=T_a_out_dsj(x);
end
%plot(t,T_g,'-..',t,T_a,':*',t,T_pv,'--+')
plot(t,T_g,'*r')
hold on
plot(t,T_a,'*y')
hold on
plot(t,T_a_out,'*g')
hold on
plot(t,T_pv,'*b')
hold on
end
xlswrite('data.xls',data)
https://jingyan.baidu.com/article/95c9d20d58ff34ec4e756197.html