高等工程热力学作业.docx

上传人:b****6 文档编号:6286166 上传时间:2023-01-05 格式:DOCX 页数:17 大小:640.98KB
下载 相关 举报
高等工程热力学作业.docx_第1页
第1页 / 共17页
高等工程热力学作业.docx_第2页
第2页 / 共17页
高等工程热力学作业.docx_第3页
第3页 / 共17页
高等工程热力学作业.docx_第4页
第4页 / 共17页
高等工程热力学作业.docx_第5页
第5页 / 共17页
点击查看更多>>
下载资源
资源描述

高等工程热力学作业.docx

《高等工程热力学作业.docx》由会员分享,可在线阅读,更多相关《高等工程热力学作业.docx(17页珍藏版)》请在冰豆网上搜索。

高等工程热力学作业.docx

高等工程热力学作业

高等工程热力学作业(编程)

第三章实际气体状态方程

第四章实际气体导出热力学性质与过程

题目:

一、用PR方程计算制冷剂R290、R600a和混合制冷剂R290/R600a:

50/50wt%的PVT性质。

二、用PR方程计算制冷剂R290、R600a和混合制冷剂R290/R600a的导出热力学性质焓和熵。

源程序:

1、牛顿迭代法求Z

functionZ=newton(A,B,Z)

err=1e-6;

forn=0:

1000

f=Z^3-(1-B)*Z^2+Z*(A-2*B-3*B^2)-(A*B-B^2-B^3);

Z=Z-f/(3*Z^2-2*(1-B)*Z+(A-2*B-3*B^2));

if(abs(f)

break

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2、求a、b、Z、v等参数函数

function[v,Z,a,b,beta]=vv(p,T)

R=8.31451;

N1=[44.096369.894.25120.1521];

N2=[58.122407.813.62900.1840];

k1=0.37464+1.54226*N1(4)-0.26992*N1(4)^2;

alpha1=(1+k1*(1-(T/N1

(2))^0.5))^2;

a1=0.45724*alpha1*R^2*N1

(2)^2/N1(3)/10^6;

aa1=0.45724*R^2*N1

(2)^2/N1(3)/10^6*2*sqrt(alpha1)*(-k1/(2*sqrt(N1

(2)*T)));

b1=0.07780*R*N1

(2)/N1(3)/10^6;

k2=0.37464+1.54226*N2(4)-0.26992*N2(4)^2;

alpha2=(1+k2*(1-(T/N2

(2))^0.5))^2;

a2=0.45724*alpha2*R^2*N2

(2)^2/N2(3)/10^6;

aa2=0.45724*R^2*N2

(2)^2/N2(3)/10^6*2*sqrt(alpha2)*(-k2/(2*sqrt(N2

(2)*T)));

b2=0.07780*R*N2

(2)/N2(3)/10^6;

a3=0.25*a1+0.5*(1-0.01)*sqrt(a1*a2)+0.25*a2;

aa3=0.25*aa1+0.5*(1-0.01)*1/2/sqrt(a1*a2)*(a1*aa2+a2*aa1)+0.25*aa2;

b3=0.5*(b1+b2);

a=[a1a2a3];

b=[b1b2b3];

beta=[aa1aa2aa3];

fori=1:

3;

A(i)=a(i)*p*10^6/(R^2*T^2);

B(i)=b(i)*p*10^6/(R*T);

Z(i)=newton(A(i),B(i),1);

vv(i)=R*T*Z(i)/p/10^6;

digits(5);

v(i)=vpa(vv(i),5);

i=i+1;

end

a=[a1a2a3];

b=[b1b2b3];

beta=[aa1aa2aa3];

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3、余函数法求ar、sr、hr

function[ar,sr,hr]=as(p,T)

[v,Z,a,b,beta]=vv(p,T);

R=8.31451;

fori=1:

3;

sr(i)=-R*log((v(i)-b(i))/v(i))+beta(i)/(2*sqrt

(2)*b(i))*log((v(i)-0.414*b(i))/(v(i)+2.414*b(i)))-R*log(v(i)/(R*T/p/10^6));

ar(i)=R*T*log((v(i)-b(i))/v(i))-a(i)/(2*sqrt

(2)*b(i))*log((v(i)-0.414*b(i))/(v(i)+2.414*b(i)))+R*T*log(v(i)/(R*T/p/10^6));

hr(i)=ar(i)+T*sr(i)+R*T*(1-Z(i));

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

4求绝对焓熵(以0℃饱和液体为标准

function[h,s]=hs(p,T)

M1=44.096;

M2=58.122;

x1=(1/M1)/(1/M1+1/M2);

x2=(1/M2)/(1/M1+1/M2);

Mm=M1*x1+M2*x2;

M=[M1M2Mm];

ps=[0.0156960.329790.47446];

T0=273.15;

R=8.31451;

c1=[-95.806.945-3.597*10^(-3)7.290*10^(-7)];

c2=[-23.916.605-3.176*10^(-3)4.981*10^(-7)];

c3=[-64.796.798-3.415*10^(-3)6.294*10^(-7)];

cps1=inline('-95.80./t+6.945-3.597*10^(-3)*t+7.290*10^(-7)*t.^2');

cps2=inline('-23.91./t+6.605-3.176*10^(-3)*t+4.981*10^(-7)*t.^2');

cps3=inline('-64.79./t+6.798-3.415*10^(-3)*t+6.294*10^(-7)*t.^2');

cph1=inline('-95.80+6.945*t-3.597*10^(-3)*t.^2+7.290*10^(-7)*t.^3');

cph2=inline('-23.91+6.605*t-3.176*10^(-3)*t.^2+4.981*10^(-7)*t.^3');

cph3=inline('-64.79+6.798*t-3.415*10^(-3)*t.^2+6.294*10^(-7)*t.^3');

Is1=quad(cps1,273.15,T)/1000;

Is2=quad(cps2,273.15,T)/1000;

Is3=quad(cps3,273.15,T)/1000;

Ih1=quad(cph1,273.15,T)/1000;

Ih2=quad(cph2,273.15,T)/1000;

Ih3=quad(cph3,273.15,T)/1000;

Is=[Is1Is2Is3];

Ih=[Ih1Ih2Ih3];

[ar,sr,hr]=as(p,T);

fori=1:

3

[ar1,sr1,hr1]=as(ps(i),T0);

ar0(i)=ar1(i);

sr0(i)=sr1(i);

hr0(i)=hr1(i);

s(i)=1*M(i)+sr0(i)+Is(i)*M(i)-R*log(p/ps(i))-sr(i);

h(i)=200*M(i)+hr0(i)+Ih(i)*M(i)-hr(i);

end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

5、主程序求v、h、s

clear

P=input('输入R600a工质压力:

P/MPa:

\n');

T=input('输入R600a工质温度:

T/K:

\n');

[v]=vv(p,T)

[h,s]=hs(p,T)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R290、R600a、R290/R600a的比体积v/(m^3/mol);

R290、R600a、R290/R600a的焓h/(J/mol);

R290、R600a、R290/R600a的熵s/(J/(mol.K);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

运行结果:

单位制:

SI

第六章气液相平衡

题目:

试用Peng-Robinson方程计算纯质R290P-T相图和溶液R290/R600a分别在p=1atm和p=10atm下的T-X相图。

纯质R290P-T相图:

源程序:

1、求纯质R290逸度系数函数

function[phi1]=phi(T1,P1,Z);

R=8.3145;

M1=44.096e-3;

Tc1=369.89;

Pc1=4.2512e6;

w1=0.1512;

Tr1=T1/Tc1;

k1=0.37464+1.54226*w1-0.26992*w1^2;

alpha1=(1+k1*(1-Tr1^0.5))^2;

a1=0.45724*alpha1*(R^2)*(Tc1^2)/Pc1;

b1=0.07780*R*Tc1/Pc1;

A1=a1*P1/((R^2)*(T1^2));

B1=b1*P1/(R*T1);

Z=newton(A1,B1,Z);

phi1=exp(Z-1-log(Z-B1)-A1*log((Z+2.414*B1)/(Z-0.414*B1))/(2*sqrt

(2)*B1));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2、主程序

p1=3e5;dp=100;

N=20000;err=1e-8;

forT1=200:

0.1:

369.89

forn=1:

N

phi1v=phi(T1,p1,1.1);

phi1L=phi(T1,p1,0.001);

ifabs(phi1v-phi1L)<=err

break

else

p1=p1+dp;

end

end

ifn==N+1

fprintf('error!

')

break;

else

plot(T1,p1/10^6,'r-');

holdon;

end

end

grid;

title('R290工质p-T相图');

xlabel('T/K');ylabel('p/MPa');

运行结果:

溶液R290/600a分别在P=1atm和10atm下的T-X相图

源程序:

1、求R290/R600a逸度系数函数

functionphimix=phimix(type,x1,T,p,Z)

R=8.3145;

x2=1-x1;

N1=[44.096369.894.25120.1521];

N2=[58.122407.813.62900.1840];

k1=0.37464+1.54226*N1(4)-0.26992*N1(4)^2;

alpha1=(1+k1*(1-(T/N1

(2))^0.5))^2;

a1=0.45724*alpha1*R^2*N1

(2)^2/N1(3)/10^6;

b1=0.07780*R*N1

(2)/N1(3)/10^6;

k2=0.37464+1.54226*N2(4)-0.26992*N2(4)^2;

alpha2=(1+k2*(1-(T/N2

(2))^0.5))^2;

a2=0.45724*alpha2*R^2*N2

(2)^2/N2(3)/10^6;

b2=0.07780*R*N2

(2)/N2(3)/10^6;

a=x1*x1*a1+2*x1*x2*(1-0.01)*sqrt(a1*a2)+x2*x2*a2;

b=x1*b1+x2*b2;

A=a*p*10^6/(R^2*T^2);

B=b*p*10^6/(R*T);

Z=newton(A,B,Z);

iftype==1

bi=b1;

sai=2*(x1*a1+x2*0.99*sqrt(a1*a2));

elseiftype==2

bi=b2;

sai=2*(x2*a2+x1*0.99*sqrt(a1*a2));

end

end

phimix=exp(bi/b*(Z-1)-log(Z-B)-A*(sai/a-bi/b)*log((Z+2.414*B)/(Z-0.414*B))/(2*sqrt

(2)*B));

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2、1atm下R290/R600a的T-x图

clear

x1=0:

0.01:

1;

x2=1-x1;

t=length(x1);

y1=x1;

y2=1-y1;

P=0.1;

n=0;

fori=1:

t

T=220;

while1

n=n+1;

fug1_l=phimix(1,x1(i),T,P,0.01);

fug1_v=phimix(1,y1(i),T,P,1.1);

fug2_l=phimix(2,x1(i),T,P,0.01);

fug2_v=phimix(2,y1(i),T,P,1.1);

k1=fug1_l/fug1_v;

k2=fug2_l/fug2_v;

y1(i)=k1*x1(i);

y2(i)=k2*x2(i);

sumy=y1(i)+y2(i);

sumy1=sumy;

ifn==1

y1(i)=k1*x1(i)/sumy;

y2(i)=k2*x2(i)/sumy;

fug1_v=phimix(1,y1(i),T,P,1.1);

fug2_v=phimix(2,y1(i),T,P,1.1);

k1=fug1_l/fug1_v;

k2=fug2_l/fug2_v;

y1(i)=k1*x1(i);

y2(i)=k2*x2(i);

sumy=y1(i)+y2(i);

end

while1

ifabs((sumy-sumy1)/sumy1)<1e-4

break

end

sumy1=sumy;

y1(i)=k1*x1(i)/sumy;

y2(i)=k2*x2(i)/sumy;

fug1_v=phimix(1,y1(i),T,P,1.1);

fug2_v=phimix(2,y1(i),T,P,1.1);

k1=fug1_l/fug1_v;

k2=fug2_l/fug2_v;

y1(i)=k1*x1(i);

y2(i)=k2*x2(i);

sumy=y1(i)+y2(i);

end

ifabs(sumy-1)<=1e-4

q(i)=T;

break

end

T=T+0.01;

end

R(:

i)=[x1(i),y1(i),T];

end

s=0;

fori=1:

t

ifR(3,i)<265

s=s+1;

L(:

s)=R(:

i);

end

end

L(:

1)=[0,0,261];

L(:

s+1)=[1,1,230.61];

plot(L(1,:

),L(3,:

),'r');holdon;

plot(L(2,:

),L(3,:

));

legend('泡点线',’露点线’);

xlabel('R290摩尔分数');

ylabel('混合工质温度/K');

title('p=1atm下,R290/R600a的T-x图');

gridon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

3、10atm下R290/R600a的T-x图

clear

x1=0:

0.01:

1;

x2=1-x1;

t=length(x1);

y1=x1;

y2=1-y1;

P=1;

n=0;

fori=1:

t

T=290;

while1

n=n+1;

fug1_l=phimix(1,x1(i),T,P,0.01);

fug1_v=phimix(1,y1(i),T,P,1.1);

fug2_l=phimix(2,x1(i),T,P,0.01);

fug2_v=phimix(2,y1(i),T,P,1.1);

k1=fug1_l/fug1_v;

k2=fug2_l/fug2_v;

y1(i)=k1*x1(i);

y2(i)=k2*x2(i);

sumy=y1(i)+y2(i);

sumy1=sumy;

ifn==1

y1(i)=k1*x1(i)/sumy;

y2(i)=k2*x2(i)/sumy;

fug1_v=phimix(1,y1(i),T,P,1.1);

fug2_v=phimix(2,y1(i),T,P,1.1);

k1=fug1_l/fug1_v;

k2=fug2_l/fug2_v;

y1(i)=k1*x1(i);

y2(i)=k2*x2(i);

sumy=y1(i)+y2(i);

end

while1

ifabs((sumy-sumy1)/sumy1)<1e-4

break

end

sumy1=sumy;

y1(i)=k1*x1(i)/sumy;

y2(i)=k2*x2(i)/sumy;

fug1_v=phimix(1,y1(i),T,P,1.1);

fug2_v=phimix(2,y1(i),T,P,1.1);

k1=fug1_l/fug1_v;

k2=fug2_l/fug2_v;

y1(i)=k1*x1(i);

y2(i)=k2*x2(i);

sumy=y1(i)+y2(i);

end

ifabs(sumy-1)<=1e-4

q(i)=T;

break

end

T=T+0.01;

end

R(:

i)=[x1(i),y1(i),T];

end

plot(R(1,:

),R(3,:

),'r');holdon;

plot(R(2,:

),R(3,:

));

legend('泡点线',’露点线’);

xlabel('R290摩尔分数');

ylabel('混合工质温度/K');

title('p=10atm下,R290/R600a的T-x图');

gridon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

运行结果:

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置:首页 > 工程科技 > 电力水利

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1