高等工程热力学作业Word文档下载推荐.docx
《高等工程热力学作业Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《高等工程热力学作业Word文档下载推荐.docx（19页珍藏版）》请在冰点文库上搜索。
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
3、余函数法求ar、sr、hr
function[ar,sr,hr]=as（p,T）
[v,Z,a,b,beta]=vv（p,T）;
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））;
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;
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）;
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）;
5、主程序求v、h、s
clear
P=input（'
输入R600a工质压力:
P/MPa:
\n'
T=input（'
输入R600a工质温度:
T/K:
[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
else
p1=p1+dp;
ifn==N+1
fprintf（'
error!
'
break;
plot（T1,p1/10^6,'
r-'
holdon;
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）
x2=1-x1;
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））;
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））;
2、1atm下R290/R600a的T-x图
x1=0:
0.01:
1;
t=length（x1）;
y1=x1;
y2=1-y1;
P=0.1;
n=0;
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;
ifabs（（sumy-sumy1）/sumy1）<
1e-4
ifabs（sumy-1）<
=1e-4
q（i）=T;
T=T+0.01;
R（:
i）=[x1（i）,y1（i）,T];
s=0;
ifR（3,i）<
265
s=s+1;
L（:
s）=R（:
i）;
L（:
1）=[0,0,261];
s+1）=[1,1,230.61];
plot（L（1,:
）,L（3,:
）,'
r'
plot（L（2,:
））;
legend（'
泡点线'
，’露点线’）;
R290摩尔分数'
混合工质温度/K'
p=1atm下，R290/R600a的T-x图'
gridon
3、10atm下R290/R600a的T-x图
P=1;
T=290;
plot（R（1,:
）,R（3,:
plot（R（2,:
p=10atm下，R290/R600a的T-x图'
gridon
升级会员 