高等工程热力学作业

2014年流体热物性部分

高等工程热力学作业

一、请用PR方程计算工质R290、R600a和混合工质R290/R600a(50/50Wt%)的

pvT性质。

PR方程的形式为：

p?RTa? v?bv(v?b)?b(v?b)a?0.45724?R2Tc2pc …… b?0.07780RTc pc??[1?k(1?Tr0.5)]2

k?0.37464?1.54226w?0.26992w2

am???xixj(1?kij)aiaj …… bm??xibi

ijik12?0.01 …………(R32/R125)

已知参数：R=8.31451;

R290： M=44.096g/mol Tc=369.89K pc=4.2512MPa w=0.1521 R600a：M=58.122g/mol Tc=407.81K pc=3.629Pa w=0.184

解题思路：

本题采用PR方程求解工质R290和R290a及其混合物的物性，已知温度和压力即可得到摩尔体积的值。其编写程序如下。 1、程序如下 (1)纯质R290

clc;clear; pc1=4.251200; tc1=369.89; R=8.314;

p=input('输入工质压力MPa)'); t=input('输入工质温度℃'); w1=0.1521; t=t+273.15; pr1=p/pc1; tr1=t/tc1;

k1=0.37464+1.54226*w1-0.26992*w1^2; al1=(1+k1*(1-tr1^0.5))^2; al=0.45724*al1*R^2*tc1^2/pc1; b1=0.07780*R*tc1/pc1; A1=0.45724*al1*pr1/tr1^2;

B1=0.07780*pr1/tr1;

z=1.0;D=1;fz=1; %给变量赋初值 eps1=0.00001;eps2=0.00001; %采用牛顿迭代法求z的值 while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B1)*z^2+z*(A1-2*B1-3*B1^2)-(A1*B1-B1^2-B1^3); fzz=3*z^2-2*(1-B1)*z+A1-2*B1-3*B1^2; y=z; z=z-fz/fzz; D=abs((z-y)/z); end

v=z*8.31451*t/p/10^6;

fprintf('该温度及压力下的工质的摩尔体积为v=%.6fm3/mol\\n',v);

(2)R600a

同R290，只是物性参数改变。 clear;clc； pc2=3.629; R=8.314; tc2=407.81;

p=input('输入工质的压力(MPa)'); t=input('输入工质的温度(℃)'); w2=0.184;

%定义计算需要的系数 t=t+273.15; pr2=p/pc2; tr2=t/tc2;

k2=0.37464+1.54226*w2-0.26992*w2^2; al2=(1+k2*(1-tr2^0.5))^2;

a2=0.45724*al2*R^2*tc2^2/pc2; b2=0.07780*R*tc2/pc2; A2=0.45724*al2*pr2/tr2^2; B2=0.07780*pr2/tr2; %定义迭代的初值 z=1.0;D=1;fz=1;

eps1=0.00001;eps2=0.00001;

%牛顿迭代法求解压缩因子Z 的方程及收敛条件 while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B2)*z^2+z*(A2-2*B2-3*B2^2)-(A2*B2-B2^2-B2^3); fzz=3*z^2-2*(1-B2)*z+A2-2*B2-3*B2^2; y=z;

z=z-fz/fzz; D=abs((z-y)/z); end

v=z*8.31451*t/p/10^6;

fprintf('该温度及压力下的工质的摩尔体积为v=%.6fm3/mol\\n',v);

(3)对于混合工质按照混合法则求解 clear;clc

p=input('输入工质的压力(MPa)'); t=input('输入工质的温度(℃)');

pc1=4.251200;tc1=369.89;w1=0.1521;t=t+273.15; pr1=p/pc1;tr1=t/tc1;

k1=0.37464+1.54226*w1-0.26992*w1^2; al1=(1+k1*(1-tr1^0.5))^2;

a1=0.45724*al1*R^2*tc1^2/pc1; b1=0.07780*R*tc1/pc1;

pc2=3.629;R=8.314;tc2=407.81;w2=0.184; t=t+273.15; pr2=p/pc2; tr2=t/tc2;

k2=0.37464+1.54226*w2-0.26992*w2^2; al2=(1+k2*(1-tr2^0.5))^2;

a2=0.45724*al2*R^2*tc2^2/pc2; b2=0.07780*R*tc2/pc2;

k12=0.01;k21=0.01;R=8.314;

am=0.3233*a1+0.1861*a2+0.4857*(1-k12)*sqrt(a1*a2); bm=0.5686*b1+0.4314*b2; Am=am*p/t^2/R^2; Bm=bm*p/t/R; %定义迭代的初值 z=1.0;D=1;fz=1;

eps1=0.00001;eps2=0.00001;

%牛顿迭代法求解压缩因子Z 的方程及收敛条件 while D>eps1&&abs(fz)>eps2

fz=z^3-(1-Bm)*z^2+z*(Am-2*Bm-3*Bm^2)-(Am*Bm-Bm^2-Bm^3); fzz=3*z^2-2*(1-Bm)*z+Am-2*Bm-3*Bm^2; y=z;

z=z-fz/fzz; D=abs((z-y)/z); end

v=z*8.31451*t/p/10^6;

fprintf('对应于该温度及压力下的该工质的比体积（m3/mol）为v=%.6fm3/mol\\n',v);

2、计算结果

按照以上编程，用PR方程分别计算工质压力为1MPa.温度为27℃时，其摩尔体积值。 （1）R290

p=1MPa.t=27℃时，摩尔体积为0.002034m/mol

3

（2）R600a

p=1MPa.t=27℃时，摩尔体积为0.001643m/mol

3（3）混合工质

p=1MPa.t=27℃时，摩尔体积为0.004596m/mol

3

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

二、请用PR方程计算工质R290、R600a和混合工质R290/R600a(50/50Wt%)的导出热力性质焓和熵。PR方程的余焓和余熵函数为：

ar?a*?a?RTlnv?bav?0.414bv?ln?RTln* vv22bv?2.414bv?b?v?0.414bv?ln?RTln* vv22bv?2.414bbr?b*?b??RTln?a ?T??理想气体比热：cp0=c0+ c1T+ c2T2 + c3T3 J/(kg.K)

R290： c0=-95.8 c1=6.945 c2=-3.597×10-3 c3=7.290×10-7 R600a： c0=-29.91 c1=6.605 c2=-3.176×10-3 c3=4.981×10-7 0oC工质的饱和压力：

R290：ps=0.15696 MPa R600a：ps=0.47446MPa R290/R600a(50/50Wt%)：ps=0.32979 MPa

解题思路：

按照国际标准以0℃饱和液体为计算基准点。

取h0℃,ps?200kJ/kg s0℃,ps?1kJ/kg?K 焓值计算

273.15hp,T?0?熵值计算

?0cdT??0pT273.15cdT?hrp,T?200?hr0℃,ps??0pT273.15c0pdT?hrp,T

273.15sp,T?0??0c0pTc0pspppdT??dT?Rln?Rln?srp,T?1?sr0℃,ps??dT?Rln?srp,T273.15T273.15TTp0pspsTc0p1、编程如下：

1）对于纯质R290 clear;clc

%输入已知参数

p=input('输入待求工质的压力(MPa)'); t=input('输入待求工质的温度(℃)'); t0=273.15; t=t+t0;

R=8.31451;

%R290工质的参数 pc=4.2512; tc=369.89;

w=0.1521; M=44.096; ps=0.15696; %理想气体比热 c0=-95.80;

c1=6.945;c2=-3.597*10^-3; c3=7.290*10^-7;

%求解余焓和余熵函数 [hr sr]=hr_s(p,t,pc,tc,w,M);

fprintf('该温度及压力下的该工质的余焓为hr=%.6f kJ/kg\\n',hr); fprintf('该温度及压力下的该工质的余熵为sr=%.6f kJ/kgK\\n',sr); [hr0 sr0]=hr_s(ps,273.15,pc,tc,w,M);

fprintf('0℃时工质饱和液体的余焓为hr=%.6f kJ/kg\\n',hr0); fprintf('0℃时工质饱和液体的余熵为sr=%.6f kJ/kgK\\n',sr0); %从273.15度到T进行积分cp0和cp0/T的值

int1=c0*t+1/2*c1*t^2+1/3*c2*t^3+1/4*c3*t^4-(c0*t0+1/2*c1*t0^2+1/3*c2*t0^3+1/4*c3*t0^4); int2=c1*t+1/2*c2*t^2+1/3*c3*t^3-c0/t^2-(c1*t0+1/2*c2*t0^2+1/3*c3*t0^3-c0/t0^2); %计算工质的焓和熵 h=200+hr0+int1/1000-hr;

s=1+sr0+int2/1000-R*log(p/ps)/M-sr; %输出计算结果

fprintf('该温度及压力下工质的焓为h=%.6f kJ/kg\\n',h); fprintf('该温度及压力下工质的熵为s=%.6f kJ/kgK\\n',s); 2) 对于纯质R600a clear;

p=input('输入待求工质的压力(MPa)'); t=input('输入待求工质的温度(℃)'); t0=273.15;t=t+t0; R=8.31451;

%R600a工质的参数 pc=3.629; tc=407.81; w=0.184; M=58.122; ps=0.47446;

%理想气体比热参数 c0=-23.91; c1=6.605;

c2=-3.176*10^-3; c3=4.981*10^-7;

%求解余焓和余熵函数 [hr, sr]=hr_s (p,t,pc,tc,w,M);

fprintf('该温度及压力下的该工质的余焓为hr=%.6f kJ/kg\\n',hr); fprintf('该温度及压力下的该工质的余熵为sr=%.6f kJ/kgK\\n',sr);

[hr0, sr0]=hr_s (ps,273.15,pc,tc,w,M);

fprintf('0℃时工质饱和液体的余焓为hr=%.6f kJ/kg\\n',hr0); fprintf('0℃时工质饱和液体的余熵为sr=%.6f kJ/kgK\\n',sr0); %从273.15度到T进行积分cp0和cp0/T的值

int1=c0*t+1/2*c1*t^2+1/3*c2*t^3+1/4*c3*t^4-(c0*t0+1/2*c1*t0^2+1/3*c2*t0^3+1/4*c3*t0^4); int2=c1*t+1/2*c2*t^2+1/3*c3*t^3-c0/t^2-(c1*t0+1/2*c2*t0^2+1/3*c3*t0^3-c0/t0^2); %计算工质的焓和熵 h=200+hr0+int1/1000-hr;

s=1+sr0+int2/1000-R*log(p/ps)/M-sr; %输出计算结果

fprintf('该温度及压力下工质的焓为h=%.6f kJ/kg\\n',h); fprintf('该温度及压力下工质的熵为s=%.6f kJ/kgK\\n',s);

3) 对于混合工质R600a/R290 W/W=50/50 clear

p=input('输入工质的压力(MPa)'); t=input('输入工质的温度(℃)'); t0=273.15; t=t+t0;

R=8.31451; ps=0.32979; M1=44.096; M2=58.122;

M=0.5686*M1+0.4314*M2; %R290理想气体比热 c01=-95.80; c11=6.945;

c21=-3.597*10^-3; c31=7.290*10^-7;

%R600a理想气体比热 c02=-23.91; c12=6.605;

c22=-3.176*10^-3; c32=4.981*10^-7;

%混合工质理想气体比热 c0=0.5686*c01+0.4314*c02; c1=0.5686*c11+0.4314*c12; c2=0.5686*c21+0.4314*c22; c3=0.5686*c31+0.4314*c32; %求解余焓和余熵函数 [hr sr]=Mix(p,t);

fprintf('该温度及压力下工质的余焓为hr=%.6f kJ/kg\\n',hr); fprintf('该温度及压力下工质的余熵为sr=%.6f kJ/kgK\\n',sr); [hr0 sr0]=Mix(ps,t0);

fprintf('0℃时工质饱和液体的余焓为hr=%.6f kJ/kg\\n',hr0); fprintf('0℃时工质饱和液体的余熵为sr=%.6f kJ/kgK\\n',sr0); %从273.15度到T进行积分cp0和cp0/T的值

int1=c0*t+1/2*c1*t^2+1/3*c2*t^3+1/4*c3*t^4-(c0*t0+1/2*c1*t0^2+1/3*c2*t0^3 +1/4*c3*t0^4);

int2=c1*t+1/2*c2*t^2+1/3*c3*t^3-c0/t^2-(c1*t0+1/2*c2*t0^2+1/3*c3*t0^3-c0/t0^2); %计算工质和焓和熵 h=200+hr0+int1/1000-hr;

s=1+sr0+int2/1000-R*log(p/ps)/M-sr; %输出计算结果

fprintf('该温度及压力下混合工质的焓为h=%.6f kJ/kg\\n',h); fprintf('该温度及压力下混合工质的熵为s=%.6f kJ/kg?K\\n',s);

%调用了一个余函数

function [hr sr]=hr_s(p,t,pc,tc,w,M) R=8.31451; tr=t/tc;

k=0.37464+1.54226*w-0.26992*w^2; al=(1+k*(1-tr^0.5))^2;

al1=-(1+k*(1-tr^0.5))*k*(1/t/tc)^0.5; a=0.45724*al*R^2*tc^2/pc*10^-6; a1=0.45724*al1*R^2*tc^2/pc*10^-6; b=0.07780*R*tc/pc*10^-6; A=a*p*10^6/(R*t)^2; B=b*p*10^6/t/R; z=1.0;D=1;fz=1;

eps1=0.00001;eps2=0.00001;

%利用牛顿迭代法求解压缩因子Z while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B)*z^2+z*(A-2*B-3*B^2)-(A*B-B^2-B^3); fzz=3*z^2-2*(1-B)*z+A-2*B-3*B^2; y=z;

z=z-fz/fzz; D=abs((z-y)/z); end

v=z*8.31451*t/p/10^6; v1=R*t/p/10^6;

hr=(t*a1-a)/(2^1.5*b)*log((v-0.414*b)/(v+2.414*b))+R*t*(1-z);

sr=a1/(2^1.5*b)*log((v-0.414*b)/(v+2.414*b))-R*log((v-b)/v)-R*log(v/v1); hr=hr/M; sr=sr/M;

以下是混合工质参数混合原则 function [hr sr]=Mix (p,t) R=8.31451;

pc1=4.2512;tc1=369.89;w1=0.1521;M1=44.096; tr1=t/tc1;

pc2=3.629;tc2=407.81;w2=0.184;M2=58.122; tr2=t/tc2;

M=0.5686*M1+0.4314*M2;

%计算工质R290PR方程的系数

k1=0.37464+1.54226*w1-0.26992*w1^2; al1=(1+k1*(1-tr1^0.5))^2;

al11=-(1+k1*(1-tr1^0.5))*k1*(1/t/tc1)^0.5; a1=0.45724*al1*R^2*tc1^2/pc1*10^-6; a11=0.45724*al11*R^2*tc1^2/pc1*10^-6; b1=0.07780*R*tc1/pc1*10^-6; %计算工质R600aPR方程的系数

k2=0.37464+1.54226*w2-0.26992*w2^2; al2=(1+k2*(1-tr2^0.5))^2;

al21=-(1+k2*(1-tr2^0.5))*k2*(1/t/tc2)^0.5; a2=0.45724*al2*R^2*tc2^2/pc2*10^-6; a21=0.45724*al21*R^2*tc2^2/pc2*10^-6; b2=0.07780*R*tc2/pc2*10^-6; %混合法则 k12=0.01;

am=0.3233*a1+0.1861*a2+0.4857*(1-k12)*sqrt(a1*a2);

am1=0.3233*a11+0.1861*a21+0.4857*(1-k12)*(1/sqrt(4*a1*a2)*(a11*a2+a1*a21)); bm=0.5686*b1+0.4314*b2; %压缩因子Z方程的系数 A=am*p*10^6/(R*t)^2; B=bm*p*10^6/t/R; z=1.0;D=1;fz=1;

eps1=0.00001;eps2=0.00001; while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B)*z^2+z*(A-2*B-3*B^2)-(A*B-B^2-B^3); fzz=3*z^2-2*(1-B)*z+A-2*B-3*B^2; y=z;

z=z-fz/fzz; D=abs((z-y)/z); end

v=z*R*t/p/10^6;

fprintf('该温度及压力下的混合制冷剂的比体积（m3/mol）为v=%.6f\\n',v); v1=R*t/p/10^6;

hr=(t*am1-am)/(2^1.5*bm)*log((v-0.414*bm)/(v+2.414*bm))+R*t*(1-z);

sr=am1/(2^1.5*bm)*log((v-0.414*bm)/(v+2.414*bm))-R*log((v-bm)/v)-R*log(v/v1); hr=hr/M; sr=sr/M;

2、计算结果如下

1）R290余熵余焓计算结果

2）R600a余熵余焓计算结果

3) 对于混合工质R600a/R290 W/W=50/50

第六章 气液相平衡

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

解题思路：

气液相平衡时，各相温度相等，各相压力相等，而且每一种组分在各相的化学势相等，每种组分在各相的逸度亦相等。由此可得

f?fil

?vi??yip??xip

则 Ki??vi?liyi?? xi??li?vi式中，fi 、?i 分别表示溶液中i组分的逸度和逸度系数，上角标v、l相应表示汽相和液相。Ki称为i组分的汽液平衡比。

i组分的逸度系数用P-R方程求得如下：

??bAln?i?i(Z?1)?ln(Z?B)?(b22B?2?xjaijjabZ?2.414B ?i)lnbZ?0.414B本题给定压力和液相组分，可假定一个温度和汽相成分。其计算框图和计算步骤如下：

开始 结束

输入p，{xi}及 打印结果 所有物性参数； 给T、{yi}赋初 是 否 ? ?调整T Kixi?1? lv计算{?}、{} ?i i 否 是 计算yi?Kixi/Kixi 计算所有的 ??yi 是否变化 lv?Ki??i/?i 计算{?iv} 是 否 是否第一次计算 迭代 yi?Kixi ?????（1）输入压力、液相成分及各组分的物性参数，并给压力及汽相成分赋初值。

（2）根据上面所提供的数据，调用计算物性而选用的状态方程中常系数的子程序，计算混合物的状态方程中的常系数。

（3）求出才常系数后，根据给定的压力及假设的温度迭代求解比容值：用液相摩尔百分数xi代入混合法则中的摩尔百分数项，迭代求出液相比容；用汽相摩尔百分数yi代入混合法则中的摩尔百分数项，迭代求出汽相比容。

（4）根据所选用的状态方程的表达式，根据（3）算出的比容，已知压力、xi及假设的p、yi，分别计算出。 （5）计算平衡比

（6）判断是否为第一次迭代。若为第一次迭代，则重新计算，并利用新的yi求方程中的系数，在迭代求解汽相比容，重复第五步。若不是第一次迭代，则坐第七步。

（7）判断是否在变化。若变化，则重新计算，类似第六步的作法，重新计算；若不变化，再判断是否等于1.若，则说米假定的温度T的初值不对，重新调整T，再返回第二步；若，则将结果输出。

程序如下

1、纯质R290p-T相图

1）求R290液相逸度系数子程序

function fl=fl_s(p,t) t=t+273.15;

pc=4.2512;tc=369.89;w=0.1521;%R290的物性参数 tr=t/tc; R=8.31451;

%求解PR方程的各个常系数 k=0.37464+1.54226*w-0.26992*w^2; al=(1+k*(1-tr^0.5))^2;

a=0.45724*al*R^2*tc^2/pc*10^-6; b=0.07780*R*tc/pc*10^-6; A=a*p*10^6/(R*t)^2; B=b*p*10^6/t/R; z=0.0001;D=1;fz=1; eps1=0.00001;eps2=0.00001; %采用牛顿迭代法求解z的值 while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B)*z^2+z*(A-2*B-3*B^2)-(A*B-B^2-B^3); fzz=3*z^2-2*(1-B)*z+A-2*B-3*B^2; y=z; z=z-fz/fzz; D=abs((z-y)/z); end

v=z*8.31451*t/p/10^6;%液相比容 v1=R*t/p/10^6;

fl=R*t*log((v-b)/v)+R*t*log(v/v1)-a/(2^1.5*b)*log((v-0.414*b)/(v+2.414*b))-(1-z);%液相逸度系数

的自然对数

2）求R290汽相逸度系数子程序

function fg=fg_s(p,t) t=t+273.15;

pc=4.2512;tc=369.89;w=0.1521;% R290的物性参数 tr=t/tc; R=8.31451;

%求解PR方程的各个常系数 k=0.37464+1.54226*w-0.26992*w^2; al=(1+k*(1-tr^0.5))^2;

a=0.45724*al*R^2*tc^2/pc*10^-6; b=0.07780*R*tc/pc*10^-6; A=a*p*10^6/(R*t)^2; B=b*p*10^6/t/R; z=1.0;D=1;fz=1;

eps1=0.00001;eps2=0.00001; %采用牛顿迭代法求解z的值 while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B)*z^2+z*(A-2*B-3*B^2)-(A*B-B^2-B^3); fzz=3*z^2-2*(1-B)*z+A-2*B-3*B^2; y=z; z=z-fz/fzz; D=abs((z-y)/z); end

v=z*8.31451*t/p/10^6; v1=R*t/p/10^6;

fg=R*t*log((v-b)/v)+R*t*log(v/v1)-a/(2^1.5*b)*log((v-0.414*b)/(v+2.414*b))-(1-z); %汽相逸度系数的自然对数

3）画R290p-t相图

clear;clf;clc; hold on

axis([-60,100,0,5]); grid p=0.5; n=0;

for t=-60:1:93 fl=1;fg=2;

while abs(fl-fg)>=0.0001 p=p+0.0001; fl=fl_s(p,t); fg=fg_s(p,t);

fprintf('计算次数n=%.f\\n',n); n=n+1; end

plot(t,p,'r.','markersize',10);hold on; end

title('R290的p-T相图'); xlabel('温度t/℃'); ylabel('压力p/MPa');

2、计算溶液R290/R600a在p=1atm下的T-x相图 1）调用汽相逸度系数子程序

function [mf1 mf2]=mf_s(p,t,y1,y2,z0) t=t+273.15; R=8.31451;

pc1=4.2512;tc1=369.89;w1=0.1521; tr1=t/tc1;

pc2=3.629;tc2=407.81;w2=0.184; tr2=t/tc2;

%计算工质R290的各个系数

k1=0.37464+1.54226*w1-0.26992*w1^2; al1=(1+k1*(1-tr1^0.5))^2; a1=0.45724*al1*R^2*tc1^2/pc1; b1=0.07780*R*tc1/pc1;

%计算工质R600a的各个系数

k2=0.37464+1.54226*w2-0.26992*w2^2; al2=(1+k2*(1-tr2^0.5))^2; a2=0.45724*al2*R^2*tc2^2/pc2; b2=0.07780*R*tc2/pc2; %混合法则 k12=0.01;

am=y1^2*a1+y2^2*a2+2*y1*y2*(1-k12)*sqrt(a1*a2); bm=y1*b1+y2*b2; %求解z

A=am*p/R^2/t^2; B=bm*p/R/t;

z=z0;D=1;fz=1; % 迭代初值，汽相z=1；液相z=0； B=bm*p/R/t; z=z0;D=1;fz=1;

eps1=0.00001;eps2=0.00001; %采用牛顿迭代法求解z的值 while D>eps1&&abs(fz)>eps2

fz=z^3-(1-B)*z^2+z*(A-2*B-3*B^2)-(A*B-B^2-B^3); fzz=3*z^2-2*(1-B)*z+A-2*B-3*B^2; y=z; z=z-fz/fzz; D=abs((z-y)/z); end

mf1=exp(b1/bm*(z-1)-log(z-B)-A/(2^1.5*B)*(2*(y1*a1+y2*(1-k12)*sqrt(a1*a2))/am-b1/bm)*log((z+2.414*B)/(z-0.414*B)));

mf2=exp(b2/bm*(z-1)-log(z-B)-A/(2^1.5*B)*(2*(y2*a2+y1*(1-k12)*sqrt(a1*a2))/am-b2/bm)*log((z+2.414*B)/(z-0.414*B)));

2）溶液R290/R600a在p=1atm下的T-x相图

clear;clf;clc;

p=0.101; %p=1atm，一个大气压下

for y1=0:0.001:1 %设定汽相求解的步长 y2=1-y1;x1=0.1;x2=1-x1;x=0;t=-60; while abs(x-1)>=0.01 t=t+0.1;

[mfl1, mfl2]=mf_s(p,t,x1,x2,0.001); %液相逸度 [mfg1, mfg2]=mf_s(p,t,y1,y2,1.001); %汽相逸度 k1=mfg1/mfl1; k2=mfg2/mfl2; x1=k1*y1/(k1*y1+k2*y2); x2=k2*y2/(k1*y1+k2*y2); x0=x; x=k1*y1+k2*y2; while abs(x-x0)>0.001

[mfl1, mfl2]=mf_s(p,t,x1,x2,0.001); k1=mfg1/mfl1; k2=mfg2/mfl2; x1=k1*y1/(k1*y1+k2*y2); x2=k2*y2/(k1*y1+k2*y2); x0=x; x=k1*y1+k2*y2; end end T=273.15+t;

fprintf(‘R290汽相组分y1=%.3f, 液相组分x1=%.3f,相平衡温度T=%.2fK\\n,y1,x1,T); axis([0,1,220,270]);

plot(x1,T,'b.','markersize',5);hold on; plot(y1,T,'r.','markersize',5);hold on; grid end

title('溶液R290/R600a在1atm下的T-x相图'); xlabel('组分x(y)'); ylabel('温度T/K');

同理，在p=10atm时，修改主程序中的压力并重新给温度赋初值，可得到同样的图像如下：