收藏
下载资源下载资源 加入VIP,免费下载

高等工程热力学作业.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