高等工程热力学作业.docx
- 文档编号:6286166
- 上传时间:2023-01-05
- 格式:DOCX
- 页数:17
- 大小:640.98KB
高等工程热力学作业.docx
《高等工程热力学作业.docx》由会员分享,可在线阅读,更多相关《高等工程热力学作业.docx(17页珍藏版)》请在冰豆网上搜索。
高等工程热力学作业
高等工程热力学作业(编程)
第三章实际气体状态方程
第四章实际气体导出热力学性质与过程
题目:
一、用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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 运行结果:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 高等 工程 热力学 作业