基于Java程序的化工热力学计算——利用状态方程预测饱和性质

f1[j]=-(z[0]-z[1])/r/t; ps[j+1]=ps[j]-f[j]/f1[j]; j++; }while((ps[j]-ps[j-1])/ps[j-1]>0.001 & j<299); System.out.println("saturated vapor pressure calculated by state equation="+String.format("%.2f", ps[j]/Math.pow(10, 6))+" MPa, ="+String.format("%.3f", ps[j]/101325)+" atm"); System.out.println("The volume of liquid is:"+String.format("%.4f", z[1]*1000)+" L/mol"); System.out.println("The volume of gas is:"+String.format("%.4f", z[0]*1000)+" L/mol"); } public static String fileToBeRead="H:/java/physicalpropertydata.xls"; public static void main(String args[]){ Scanner sz=new Scanner(System.in); double tc,pc,w,t,x,p=0; double[] h=new double[5]; srkps equation=new srkps(); tc=305.4;pc=48.8;w=0.099; pc=pc*100000; System.out.println("tc="+tc+" K, pc="+pc/Math.pow(10,6)+" MPa, w="+w); System.out.println("please input the temperature of the system(K)"); t=sz.nextDouble(); h[0]=-6.34307;h[1]=1.0163;h[2]=-1.19116;h[3]=-2.03539;
0 23.708 23.812 0.0844 0.6508
10 29.930 30.074 0.0928 0.4849
30 45.974 45.935 0.1369 0.2290
x[0]=x[1]=x[2]=-b/(3*a); else {
if(L>0) {
y1=A*b+(3*a*(-B+Math.sqrt(L)))/2; y2=A*b+(3*a*(-B- Math.sqrt(L)))/2; if(y1>0)
gs1= Math.pow(y1,1.0/3); else
gs1=- Math.pow(Math.abs(y1),1.0/3); if(y2>0)
x=1-t/tc; p=Math.exp(1/(1-x)*(h[0]*x+h[1]* Math.pow(x,1.5)+h[2]* Math.pow(x,3)+h[3]* Math.pow(x,6))+ Math.log(pc)); System.out.println("saturated vapor pressure calculated by antoine equation="+String.format("%.2f", p/Math.pow(10,6))+" MPa");
do{ e[0]=1.0;e[1]=-(r*t/ps[j]);e[2]=a/ps[j]-b*r*t/ps[j]-
Math.pow(b,2);e[3]=-a*b/ps[j]; z=root.equa(e,z);
f[j]=Math.log((z[0]-b)/(z[1]-b))-a/(b*r*t)*Math.log(z[0]*(z[1]+b)/z[1]/ (z[0]+b))-ps[j]*(z[0]-z[1])/r/t;
sz.close(); equation.srkeos(tc,pc,w,t,p); } }
结果：
实测压力 /atm
计算压力 /atm
饱和液体体积 /(L/mol)
饱和气体体积 /(L/mol)
-40 7.748 7.696 0.0669 2.1531
-20 14.154 14.153 0.0737 1.1615
gs2= Math.pow(y2,1.0/3); else
gs2=- Math.pow(Math.abs(y2),1.0/3); p1=(-2*b+gs1+gs2)/(6*a); p2= Math.sqrt(3)*(gs1-gs2)/(6*a);
x[0]=(-b-gs1-gs2)/(3*a); x[1]=p1+p2*Math.sqrt(-1.0); x[2]=p1-p2*Math.sqrt(-1.0); } else if(L==0) { k=B/A; x[0]=-b/a+k; x[1]=x[2]=-k/2; } else { t=(2*A*b-3*a*B)/(2* Math.pow(A,1.5)); rad= Math.acos(t); x[0]=(-b-2*Math.sqrt(A)* Math.cos(rad/3))/(3*a); x[1]=(-b+ Math.sqrt(A)*( Math.cos(rad/3)+ Math.sqrt(3)* Math.sin(rad/3)))/(3*a); x[2]=(-b+ Math.sqrt(A)*( Math.cos(rad/3)- Math.sqrt(3)* Math.sin(rad/3)))/(3*a); } } Arrays.sort(x); for(i=0;i<3;i++) if(x[i]>0) {r[j]=x[i];j++;} if(j==0) z[0]=r[0]; else z[0]=r[j-1]; z[1]=r[0]; return z; } void srkeos(double tc, double pc, double w,double t,double p) { double a,b,r=8.3146,m; double[] e=new double[4],z=new double[2],ps=new double[300],f=new double[300],f1=new double[300]; int j=0; sat root=new sat(); ps[0]=p; m=0.48+1.574*w-0.176* Math.pow(w,2.0); a=0.42747* Math.pow(r,2.0)* Math.pow(tc,2.0)/pc* Math.pow((1+m*(1Math.pow(t/tc,0.5))),2.0); b=0.08664*r*tc/pc;
例题参考书：《化工热力学》，陈光进等。 【例 2-5】用 SRK 方程计算某温度下乙烯的饱和蒸汽压和液体及饱和蒸汽的摩尔 体积。 解： 按照框图编写程序。
【习题-8】利用附录Ⅵ中的立方方程根的求解子程序，选用适当的状态方程，计 算乙烷在-40℃、-20℃、0℃、10℃、30℃时的饱和蒸汽压和液体及饱和蒸汽的 摩尔体积。已知乙烷在-40℃、-20℃、0℃、10℃、30℃时的饱和蒸汽压实测值 分别为7.748atm、14.154 atm、23.708 atm、29.930 atm、45.974 atm。
程序：
import java.util.Arrays; import java.util.Scanner; public class srkps {
double[] equa(double[] f,,b=f[1],c=f[2],d=f[3]; double A,B,C,L,k,t,rad; double gs1,gs2,y1,y2,p1,p2; double[] x=new double[3],r=new double[4]; int i,j=0; r[0]=0; A=b*b-3*a*c; B=b*c-9*a*d; C=c*c-3*b*d; L=B*B-4*A*C; if(A==0&&B==0)