马特拉并
2009-05-18, 10:41
主要问题是画图出错，不知道错在哪里，代码如下：
function f = therm_matrix
%therm_matrix.m
'凝汽式发电厂原则性热力系统计算'
'第一部分：输入结构参数'
[n,m,rh,no_ml,h_style,beox_style,cooler,l,scat] =structuredata
'第二部分：输入设备参数'
[Pe,p0,t0,pc,hc,pmin,tmin]=turbinedata;
[pexh,texh,ploss,dt_water,dt_scatter,pout_pu,eff_pu,t_fw]=reheaterdata;
[D_b,p_b,t_b,h_b,eff_b,prh1,trh1,prh2,trh2,eff_dprh,c_blow,c_lost]=boilerdata;
[c1,c2,c3,a_xj,p_xj,h_xj,a_sg1,h_sg1,a_sg2,h_sg2,a_sg3,h_sg3,p_sg,eff,p_ma,t_ma,pout_cp,eff_m,eff_g,H_deox,a_0,p_drum,p_bl,x_bl,h_cp,]=otherdata;
'第三部分：计算'
'1绘制蒸汽膨胀线'
' -----计算抽汽比焓-----'
for j= 1:n
if j < n
[vexh(j),hexh(j),sexh(j)] =pt(pexh(j),texh(j));
else
hexh(n)=texh(n);[texh(n),vexh(n),sexh(n),xexh(n)]=ph(pexh(n),hexh(n));
end
end
'-----计算各位置蒸汽状态-----'
[v0,h0,s0]=ptsteam(p0,t0);
p0p=(1-c1)*p0;
h0p=h0;
[t0p,v0p,s0p,x0p]=ph(p0p,h0p);
[vmin,hmin,smin]=ptsteam(pmin,tmin);
pminp=pmin*(1-c2);
hminp=hmin;
[tminp,vminp,sminp,xminp]=ph(pminp,hminp);
hlin=hexh(no_ml);
plin=(1-c3)*pexh(no_ml);
[tlin,vlin,slin,xlin]=ph(plin,hlin);
[tc,vc,sc,xc,]=ph(pc,hc);
'----绘制 H--S 图----'
p(1)=p0;
p(2)=p0p;
t(1)=t0;
t(2)=t0p;
s(1)=s0;
s(2)=s0p;
h(1)=h0;
h(2)=h0p;
for j = 3:rh +2
p(j)=pexh(j-2);
t(j)=texh(j-2);
s(j)=sexh(j-2);
h(j)=hexh(j-2);
end
p(rh+3)=pmin;
p(rh+4)=pminp;
t(rh+3)=tmin;
t(rh+4)=tminp;
s(rh+3)=smin;
s(rh+4)=sminp;
h(rh+3)=hmin;
h(rh+4)=hminp;
for j= rh+5:no_ml+4
p(j)=pexh(j-4);
t(j)=texh(j-4);
s(j)=sexh(j-4);
h(j)=hexh(j-4);
end
p(no_ml+5)=plin;
t(no_ml+5)=tlin;
s(no_ml+5)=slin;
h(no_ml+5)=hlin;
for j=no_ml+6:n+4
p(j)=pexh(j-5);
t(j)=texh(j-5);
s(j)=sexh(j-5);
h(j)=hexh(j-5);
end
p(n+5)=pexh(n);
t(n+5)=texh(n);
s(n+5)=sexh(n);
h(n+5)=hexh(n);
x(n+5)=xexh(n);
p(n+6)=pc;
t(n+6)=tc;
s(n+6)=sc;
h(n+6)=hc;
x(n+6)=xc;
p
t
h
s
x
for j=1:(rh+2)
shigh(j)=s(j);hhigh(j)=h(j);
end
plot(shigh,hhigh,'ro')
axis([6 8 2000 4000])
grid on
hold on
plot(shigh,hhigh,'-')
axis([6 8 2000 4000])
xlabel('s[kJ/(kg.K)]')
ylabel('h[kJ/kg]')
hold on
for j =1:n-rh+4
slow(j)=s(j+(rh+2));hlow(j)=h(j+(rh+2));
end
plot(slow,hlow,'ro');axis([6 8 2000 4000]);grid on ;hold on ;
plot(slow,hlow,'-');axis([6 8 2000 4000]); hold on ;
for j=1:n+6
p1=p(j);p2=p(j);h1=h(j)-30;h2=h(j)+30;
[t1,v1,s1,x1]=ph(p1,h1);[t2,v2,s2,x2]=ph(p2,h2);
sline=[s1 s2];hline=[h1,h2];
plot(sline,hline,'-');axis([6 8 2000 4000]);hold on ;
x=j~=2;y=j~=rh+4;z=j~=no_ml+5;L=x&y&z;
if L>0
str=num2str(p(j));str=strcat('p',str);text(s2-0.05,h2+50,str);
else
str=num2str(p(j));str=strcat('p',str);text(s2,h2,str);
end
str=num2str(xexh(n));str=strcat('x=',str);text((sexh(n)+0.02),hexh(n)-10,str);
str=num2str(xc);str=strcat('x=',str);text((sc+0.02),hc-10,str);
end
title('机组蒸汽膨胀过程线');
'2汇总汽水参数表'
'-----计算各加热器压力，饱和温度，饱和水比焓，加热器出口水温，加热器出口水比焓-----'
for j=1:n
p_heater(j)=pexh(j)*(1-ploss(j)/100);
ts_heater(j)=tsaturation(p_heater(j));
[vs_heater(j),hs_heater(j),ss_heater(j)]=ptwater(p_heater(j),ts_heater(j));
t_water(j)=ts_heater(j)-dt_water(j);
if j<m
p_water(j)=pout_pu;
elseif j>m
p_water(j)=pout_cp;
else
p_water(j)=p_heater(j);
end
[v_water(j),h_water(j),s_water(j)]=ptwater(p_water(j),t_water(j));
end
h_scatter(m)=h_water(m);
h_water(n+1)=h_cp;
'-----计算给水泵焓升-----'
pin_pu=H_deox*9.8/1000*10+p_heater(m)*0.94;
[vin_pu,hin_pu,sin_pu]=ptwater(pin_pu,ts_heater(m));
[tout_pu,vout_ou,hout_pu,xout_pu]=ps(pout_pu,sin_pu);
dh_pu=(hout_pu-hin_pu)/eff_pu;
'------计算疏水冷却器疏水温度疏水比焓------'
for j=1:n
if h_style(j) == 0
if cooler(j)==1
if beox_style(j)==1
hp_water(j+1)=h_water(j+1)+dh_pu;
[tp_water(j+1),vp_water(j+1),sp_water(j+1),xp_water(j+1)]=ph(pout_pu,hp_water(j+1));
t_scatter(j)=tp_water(j+1)+dt_scatter(j);
else
t_scatter(j)=t_water(j+1)+dt_scatter(j);
end
[v_scatter(j),h_scatter(j),s_scatter(j)]=ptwater(p_heater(j),t_scatter(j));
else
h_scatter(j)=hs_heater(j);
end
else
end
end
%输出数据
pexh
texh
hexh
sexh
ploss
p_heater
ts_heater
hs_heater
ss_heater
dt_water
t_water
h_water
s_water
h_scatter
'3.回热抽汽量计算'
'---锅炉排污利用系统计算---'
a_0p=a_0+a_sg1+a_sg2+a_sg3;a_b=a_0p/(1-c_lost);a_l=c_lost*a_b;a_bl=c_blow*a_b;
a_fw=a_b+a_bl;
ts_bp=tsaturation(p_drum);[v_bp,h_bp,s_bp]=ptwater(p_drum,ts_bp);[t_f,v_f,h_f,s_f]=px(p_bl,x_bl);
ts_bl=tsaturation(p_bl);[v_blh,h_blp,s_blp]=ptwater(p_bl,ts_bl);
a_f=a_bl*(eff*h_bp-h_blp)/(h_f-h_blp);a_blp=a_bl-a_f;a_ma=a_l+a_blp;
[v_ma,h_ma,s_ma]=ptwater(p_ma,t_ma);
[v_blpp,h_blpp,s_blpp]=ptwater(p_bl,t_ma+8);
h_map=a_blp*x_bl/a_ma*(h_blp-h_blpp)+h_ma;
ts_sg=tsaturation(p_sg);
[v_sg3p,h_sg3p,s_sg3p]=ptwater(p_sg,ts_sg);
Da=[0 0;0 0;0 0;a_f a_sg1;0 0;0 0;a_sg2 0;a_sg3 0];
h_Da=[0 0;0 0;0 0;h_f h_sg1;0 0;0 0;h_sg2 0;h_sg3 0];
'----计算各项抽汽系数----'
for j=1:n
if h_style(j)
'汇集式加热器';
q(j)=hexh(j)-h_water(j+1);
r(j)=h_scatter(j-1)-h_water(j+1);
tao(j)=h_water(j)-h_water(j+1);
for k=1:l(j);
q_Da(j,k)=h_Da(j,k)-h_water(j+1);
end
else
'表面式加热'
if j ==1
q(j)=hexh(j)-h_scatter(j);
r(j)=0;
tao(j)=h_water(j)-h_water(j+1);
else
q(j)=hexh(j)-h_scatter(j);
if h_style(j-1)
r(j)=0;
else
r(j)=h_scatter(j-1)-h_scatter(j);
end
tao(j)=h_water(j)-h_water(j+1);
end
for k=1:l(j)
if j==n
q_Da(j,k)=h_Da(j,k)-h_sg3p;
else
q_Da(j,k)=h_Da(j,k)-h_scatter(j);
end
end
end
end
q
r
tao
for j=1:n
for k =1:n
if j>k
A(j,k)=scat(j,k)*eff*r(j)+(not(scat(j,k)))*tao(j)+scat(j,k)*h_style(j)*(1-eff)*tao(j);
elseif j == k
A(j,k)=eff*q(j)+h_style(j)*(1-eff)*tao(j);
else
A(j,k)=0;
end
end
end
A
for j=1:n
EDa=0
for k=1:j-1
for i = 1:l(k)
EDa=EDa+Da(k,i);
end
end
T(j)=(a_fw-EDa)*(tao(j)-beox_style(j)*dh_pu);
end
T=T'
for j=1:n
B(j)=0;
for k=1:l(j)
B(j)=B(j)+Da(j,k)*(eff*q_Da(j,k)+h_style(j)*(1-eff)*tao(j));
end
end
B=B'
aexh= A\(T-B)
aexh(4)=aexh(4)+a_xj;
'---汽轮机通流部分流量平衡---'
a_sum=0;
for j= 1:n
a_sum=a_sum+aexh(j);
end
a_c=1-a_sum;
%用凝汽器质量平衡验算
a_cp=a_fw-a_sum-a_sg1-a_sg2-a_sg3-a_ma-a_f;
if abs(a_c-a_cp)<0.01
'flux equal'
else
'flux unequal'
end
'4.汽轮机汽耗量计算'
dq_rh=hmin-hexh(rh);
dh_rh=h0-hc+dq_rh;
for j=1:rh
w(j)=aexh(j)*(h0-hexh(j));
end
for j=rh+1:n
w(j)=aexh(j)*(h0-hexh(j)+dq_rh);
end
W=0;
for j=1:n
W=W+w(j);
end
W=W+a_c*dh_rh;D0=3600*Pe/eff_m/eff_g/W
'5热经济指标计算'
'---汽轮机发电机组热经济指标---'
a_rh=1;
for j=1:rh
a_rh=a_rh-aexh(j);
end
h_fw=h_water(1);
Q0=a_0p*D0*(h0-h_fw)+a_rh+a_f*D0*(h_f-h_fw)+a_ma*D0*(h_map-h_fw)
q0=Q0/Pe
eff_el=3600/q0
'---发电厂热经济指标---'
[vrh2,hrh2,srh2]=ptsteam(prh2,trh2);
[vrh1,hrh1,srh1]=ptsteam(prh1,trh1);
dq_rhp=hrh2-hrh1;
Qb=a_b*D0*(h_b-h_fw)+a_rh*D0*dq_rhp+a_bl*D0*(h_bp-h_fw)
eff_p=Q0/Qb
eff_cp=eff_b*eff_p*eff_el
q_cp=3600/eff_cp
bs=0.123/eff_cp
谢谢帮忙
function f = therm_matrix
%therm_matrix.m
'凝汽式发电厂原则性热力系统计算'
'第一部分：输入结构参数'
[n,m,rh,no_ml,h_style,beox_style,cooler,l,scat] =structuredata
'第二部分：输入设备参数'
[Pe,p0,t0,pc,hc,pmin,tmin]=turbinedata;
[pexh,texh,ploss,dt_water,dt_scatter,pout_pu,eff_pu,t_fw]=reheaterdata;
[D_b,p_b,t_b,h_b,eff_b,prh1,trh1,prh2,trh2,eff_dprh,c_blow,c_lost]=boilerdata;
[c1,c2,c3,a_xj,p_xj,h_xj,a_sg1,h_sg1,a_sg2,h_sg2,a_sg3,h_sg3,p_sg,eff,p_ma,t_ma,pout_cp,eff_m,eff_g,H_deox,a_0,p_drum,p_bl,x_bl,h_cp,]=otherdata;
'第三部分：计算'
'1绘制蒸汽膨胀线'
' -----计算抽汽比焓-----'
for j= 1:n
if j < n
[vexh(j),hexh(j),sexh(j)] =pt(pexh(j),texh(j));
else
hexh(n)=texh(n);[texh(n),vexh(n),sexh(n),xexh(n)]=ph(pexh(n),hexh(n));
end
end
'-----计算各位置蒸汽状态-----'
[v0,h0,s0]=ptsteam(p0,t0);
p0p=(1-c1)*p0;
h0p=h0;
[t0p,v0p,s0p,x0p]=ph(p0p,h0p);
[vmin,hmin,smin]=ptsteam(pmin,tmin);
pminp=pmin*(1-c2);
hminp=hmin;
[tminp,vminp,sminp,xminp]=ph(pminp,hminp);
hlin=hexh(no_ml);
plin=(1-c3)*pexh(no_ml);
[tlin,vlin,slin,xlin]=ph(plin,hlin);
[tc,vc,sc,xc,]=ph(pc,hc);
'----绘制 H--S 图----'
p(1)=p0;
p(2)=p0p;
t(1)=t0;
t(2)=t0p;
s(1)=s0;
s(2)=s0p;
h(1)=h0;
h(2)=h0p;
for j = 3:rh +2
p(j)=pexh(j-2);
t(j)=texh(j-2);
s(j)=sexh(j-2);
h(j)=hexh(j-2);
end
p(rh+3)=pmin;
p(rh+4)=pminp;
t(rh+3)=tmin;
t(rh+4)=tminp;
s(rh+3)=smin;
s(rh+4)=sminp;
h(rh+3)=hmin;
h(rh+4)=hminp;
for j= rh+5:no_ml+4
p(j)=pexh(j-4);
t(j)=texh(j-4);
s(j)=sexh(j-4);
h(j)=hexh(j-4);
end
p(no_ml+5)=plin;
t(no_ml+5)=tlin;
s(no_ml+5)=slin;
h(no_ml+5)=hlin;
for j=no_ml+6:n+4
p(j)=pexh(j-5);
t(j)=texh(j-5);
s(j)=sexh(j-5);
h(j)=hexh(j-5);
end
p(n+5)=pexh(n);
t(n+5)=texh(n);
s(n+5)=sexh(n);
h(n+5)=hexh(n);
x(n+5)=xexh(n);
p(n+6)=pc;
t(n+6)=tc;
s(n+6)=sc;
h(n+6)=hc;
x(n+6)=xc;
p
t
h
s
x
for j=1:(rh+2)
shigh(j)=s(j);hhigh(j)=h(j);
end
plot(shigh,hhigh,'ro')
axis([6 8 2000 4000])
grid on
hold on
plot(shigh,hhigh,'-')
axis([6 8 2000 4000])
xlabel('s[kJ/(kg.K)]')
ylabel('h[kJ/kg]')
hold on
for j =1:n-rh+4
slow(j)=s(j+(rh+2));hlow(j)=h(j+(rh+2));
end
plot(slow,hlow,'ro');axis([6 8 2000 4000]);grid on ;hold on ;
plot(slow,hlow,'-');axis([6 8 2000 4000]); hold on ;
for j=1:n+6
p1=p(j);p2=p(j);h1=h(j)-30;h2=h(j)+30;
[t1,v1,s1,x1]=ph(p1,h1);[t2,v2,s2,x2]=ph(p2,h2);
sline=[s1 s2];hline=[h1,h2];
plot(sline,hline,'-');axis([6 8 2000 4000]);hold on ;
x=j~=2;y=j~=rh+4;z=j~=no_ml+5;L=x&y&z;
if L>0
str=num2str(p(j));str=strcat('p',str);text(s2-0.05,h2+50,str);
else
str=num2str(p(j));str=strcat('p',str);text(s2,h2,str);
end
str=num2str(xexh(n));str=strcat('x=',str);text((sexh(n)+0.02),hexh(n)-10,str);
str=num2str(xc);str=strcat('x=',str);text((sc+0.02),hc-10,str);
end
title('机组蒸汽膨胀过程线');
'2汇总汽水参数表'
'-----计算各加热器压力，饱和温度，饱和水比焓，加热器出口水温，加热器出口水比焓-----'
for j=1:n
p_heater(j)=pexh(j)*(1-ploss(j)/100);
ts_heater(j)=tsaturation(p_heater(j));
[vs_heater(j),hs_heater(j),ss_heater(j)]=ptwater(p_heater(j),ts_heater(j));
t_water(j)=ts_heater(j)-dt_water(j);
if j<m
p_water(j)=pout_pu;
elseif j>m
p_water(j)=pout_cp;
else
p_water(j)=p_heater(j);
end
[v_water(j),h_water(j),s_water(j)]=ptwater(p_water(j),t_water(j));
end
h_scatter(m)=h_water(m);
h_water(n+1)=h_cp;
'-----计算给水泵焓升-----'
pin_pu=H_deox*9.8/1000*10+p_heater(m)*0.94;
[vin_pu,hin_pu,sin_pu]=ptwater(pin_pu,ts_heater(m));
[tout_pu,vout_ou,hout_pu,xout_pu]=ps(pout_pu,sin_pu);
dh_pu=(hout_pu-hin_pu)/eff_pu;
'------计算疏水冷却器疏水温度疏水比焓------'
for j=1:n
if h_style(j) == 0
if cooler(j)==1
if beox_style(j)==1
hp_water(j+1)=h_water(j+1)+dh_pu;
[tp_water(j+1),vp_water(j+1),sp_water(j+1),xp_water(j+1)]=ph(pout_pu,hp_water(j+1));
t_scatter(j)=tp_water(j+1)+dt_scatter(j);
else
t_scatter(j)=t_water(j+1)+dt_scatter(j);
end
[v_scatter(j),h_scatter(j),s_scatter(j)]=ptwater(p_heater(j),t_scatter(j));
else
h_scatter(j)=hs_heater(j);
end
else
end
end
%输出数据
pexh
texh
hexh
sexh
ploss
p_heater
ts_heater
hs_heater
ss_heater
dt_water
t_water
h_water
s_water
h_scatter
'3.回热抽汽量计算'
'---锅炉排污利用系统计算---'
a_0p=a_0+a_sg1+a_sg2+a_sg3;a_b=a_0p/(1-c_lost);a_l=c_lost*a_b;a_bl=c_blow*a_b;
a_fw=a_b+a_bl;
ts_bp=tsaturation(p_drum);[v_bp,h_bp,s_bp]=ptwater(p_drum,ts_bp);[t_f,v_f,h_f,s_f]=px(p_bl,x_bl);
ts_bl=tsaturation(p_bl);[v_blh,h_blp,s_blp]=ptwater(p_bl,ts_bl);
a_f=a_bl*(eff*h_bp-h_blp)/(h_f-h_blp);a_blp=a_bl-a_f;a_ma=a_l+a_blp;
[v_ma,h_ma,s_ma]=ptwater(p_ma,t_ma);
[v_blpp,h_blpp,s_blpp]=ptwater(p_bl,t_ma+8);
h_map=a_blp*x_bl/a_ma*(h_blp-h_blpp)+h_ma;
ts_sg=tsaturation(p_sg);
[v_sg3p,h_sg3p,s_sg3p]=ptwater(p_sg,ts_sg);
Da=[0 0;0 0;0 0;a_f a_sg1;0 0;0 0;a_sg2 0;a_sg3 0];
h_Da=[0 0;0 0;0 0;h_f h_sg1;0 0;0 0;h_sg2 0;h_sg3 0];
'----计算各项抽汽系数----'
for j=1:n
if h_style(j)
'汇集式加热器';
q(j)=hexh(j)-h_water(j+1);
r(j)=h_scatter(j-1)-h_water(j+1);
tao(j)=h_water(j)-h_water(j+1);
for k=1:l(j);
q_Da(j,k)=h_Da(j,k)-h_water(j+1);
end
else
'表面式加热'
if j ==1
q(j)=hexh(j)-h_scatter(j);
r(j)=0;
tao(j)=h_water(j)-h_water(j+1);
else
q(j)=hexh(j)-h_scatter(j);
if h_style(j-1)
r(j)=0;
else
r(j)=h_scatter(j-1)-h_scatter(j);
end
tao(j)=h_water(j)-h_water(j+1);
end
for k=1:l(j)
if j==n
q_Da(j,k)=h_Da(j,k)-h_sg3p;
else
q_Da(j,k)=h_Da(j,k)-h_scatter(j);
end
end
end
end
q
r
tao
for j=1:n
for k =1:n
if j>k
A(j,k)=scat(j,k)*eff*r(j)+(not(scat(j,k)))*tao(j)+scat(j,k)*h_style(j)*(1-eff)*tao(j);
elseif j == k
A(j,k)=eff*q(j)+h_style(j)*(1-eff)*tao(j);
else
A(j,k)=0;
end
end
end
A
for j=1:n
EDa=0
for k=1:j-1
for i = 1:l(k)
EDa=EDa+Da(k,i);
end
end
T(j)=(a_fw-EDa)*(tao(j)-beox_style(j)*dh_pu);
end
T=T'
for j=1:n
B(j)=0;
for k=1:l(j)
B(j)=B(j)+Da(j,k)*(eff*q_Da(j,k)+h_style(j)*(1-eff)*tao(j));
end
end
B=B'
aexh= A\(T-B)
aexh(4)=aexh(4)+a_xj;
'---汽轮机通流部分流量平衡---'
a_sum=0;
for j= 1:n
a_sum=a_sum+aexh(j);
end
a_c=1-a_sum;
%用凝汽器质量平衡验算
a_cp=a_fw-a_sum-a_sg1-a_sg2-a_sg3-a_ma-a_f;
if abs(a_c-a_cp)<0.01
'flux equal'
else
'flux unequal'
end
'4.汽轮机汽耗量计算'
dq_rh=hmin-hexh(rh);
dh_rh=h0-hc+dq_rh;
for j=1:rh
w(j)=aexh(j)*(h0-hexh(j));
end
for j=rh+1:n
w(j)=aexh(j)*(h0-hexh(j)+dq_rh);
end
W=0;
for j=1:n
W=W+w(j);
end
W=W+a_c*dh_rh;D0=3600*Pe/eff_m/eff_g/W
'5热经济指标计算'
'---汽轮机发电机组热经济指标---'
a_rh=1;
for j=1:rh
a_rh=a_rh-aexh(j);
end
h_fw=h_water(1);
Q0=a_0p*D0*(h0-h_fw)+a_rh+a_f*D0*(h_f-h_fw)+a_ma*D0*(h_map-h_fw)
q0=Q0/Pe
eff_el=3600/q0
'---发电厂热经济指标---'
[vrh2,hrh2,srh2]=ptsteam(prh2,trh2);
[vrh1,hrh1,srh1]=ptsteam(prh1,trh1);
dq_rhp=hrh2-hrh1;
Qb=a_b*D0*(h_b-h_fw)+a_rh*D0*dq_rhp+a_bl*D0*(h_bp-h_fw)
eff_p=Q0/Qb
eff_cp=eff_b*eff_p*eff_el
q_cp=3600/eff_cp
bs=0.123/eff_cp
谢谢帮忙