C 此程序是运用peng-robinson方程计算气体的逸度系数faita(1),faita(2)，是利用牛顿迭代法求解
C 然后通过origin拟合得到这些数据N1,N2,m1,m2,K1,K2
C 求解得到f1,f2,x1的过程,首先利用牛顿迭代法求解IAST理论得到x1, 然后利用表达式求解f1,f2
PROGRAM MAIN
LOGICAL converged
integer S N SQ
parameter(EPS=1.0E-06,R=83.14,epsilon=1-2**(0.5),sigma=1+2**(0.5)
*)
dimension faita(7),ns(4),ma(4),K(4),M(7),omiga(7),
*Tc(7),Pc(7),b(7),beta(7),alpha(7),q(7),as(7),Z1(7)
real alpha,beta,Tr,omiga,M,rho(7),q,a,b,SV,I(7),ZM,SSA,SSB,SSC,
*SSD,T,f1,f2
DATA M/16.043,30.070,44.097, 58.123,72.150, 44.01, 28.01/
DATA omiga/0.012,0.100,0.152,0.2,0.252,0.224,0.048/
DATA Tc/190.6,305.3,369.8,425.1,469.7,304.2,132.9/
DATA Pc/45.99,48.72,42.48,37.96,33.70,73.83,34.99/
READ(*,*)P1
C 计算peng-robinson方程中需要的一些参量
DO 3 N=1,7
P=P1/100
SV=(1-(T/Tc(N))**(0.5))**2
alpha(N)=(1+(0.37464+1.54226*omiga(N)-0.26992*omiga(N)**2)*SV)
as(N)=((0.45724*R**2*Tc(N)**2)/Pc(N))*alpha(N)
b(N)=0.07779*R*Tc(N)/Pc(N)
beta(N)=b(N)*P/(R*T)
q(N)=as(N)*alpha(N)*P/(R**2*T**2)
3 CONTINUE
C 运用牛顿迭代法开始求解压缩因子的方程
T=298.0
IF(N.LT.7.AND.N.GT.1)THEN
ZM=1
DO WHILE(.not.converged.and.ABS(MSD(ZM,beta(N),q(N))).GT.0)
SSA=MS(ZM,beta(N),q(N))
SSB=MSD(ZM,beta(N),q(N))
ZM=ZM-SSA/SSB
converged=(abs(SSA)/(SSB)).LT.EPS
END DO
ELSE
END IF
do 4 N=1,7
Z1(N)=ZM
4 CONTINUE
C 求解密度和逸度
DO 5 N=1,7
I(N)=(1/(sigma-epsilon))*log((Z1(N)+sigma*beta(N))/
*(Z1(N)+epsilon*beta(N)))
rho(N)=(P/(Z1(N)*R*T))*M(N)*EPS
faita(N)=SVF(Z1(N),beta,q(N),I(N))
5 CONTINUE
C 运用IAST理论以及上面得到的逸度系数，开始求解IAST的方程
S=0
X=1
SSC=iast(X,P,y1,faita(1),faita(2))
SSD=diast(X,P,y1,faita(1),faita(2))
DO WHILE(.not.converged.and.ABS(SSD).GT.0)
X=X-SSC/SSD
S=S+1
converged=(abs(SSC/SSD).LT.EPS)
END DO
x1=X
f1=P*y1*faita(1)/x1
f2=P*(1-y1)*faita(2)/(1-x1)
20 format(1x,3F9.4)
write(20,*)f1,f2,x1
end
function MS(SA,SB,SC)
MS=SA**3-(1-SB)*SA**2+(SC-2*SB-3*SB
***2)*SA-(SC*SB-SB**2-SB**3)
end
function MSD(SD,SE,SF)
MSD=3*SD**2-2*(1-SE)*SD+(SF-2*SE-
*3*SE**2)
end
function iast(MADD,P,YY,SG1,SG2)
dimension ns(4),ma(4),K(4)
REAL MADD
DATA ns/2.0,3.0,4.0,5.0/
DATA ma/1.0,2.0,3.0,4.0/
DATA K/5.0,6.0,7.0,8.0/
f1=P*YY*SG1/MADD
f2=P*(1-YY)*SG2/(1-MADD)
iast=ns(1)*LOG(1+K(1)*f1**ma(1))+ns(2)*LOG
*(1+K(2)*f1**ma(2))-ns(3)*LOG(1+K(3)*f2**ma(3))+ns(4)*LOG(1+
*K(4)*f2**ma(4))
end
function diast(MAD,S,Y,SG3,SG4)
dimension ns(4),ma(4),K(4)
REAL MAD
DATA ns/2.0,3.0,4.0,5.0/
data ma/1.0,2.0,3.0,4.0/
data K/5.0,6.0,7.0,8.0/
f1=S*Y*SG3/MAD
f2=S*(1-Y)*SG4/(1-MAD)
diast=(ma(1)*ns(1)*K(1)*f1*(ma(1)-1))/(1+
*K(1)*f1**ma(1))+(ma(2)*ns(2)*K(2)*f1*(ma(2)-1))/(1+K(2)*f1**ma(2))
*-(ma(3)*ns(3)*K(3)*f2*(ma(3)-1))/(1+K(3)*f2**ma(3))+
*(ma(4)*ns(4)*K(4)*f2*(ma(4)-1))/(1+K(4)*f2**ma(4))
END
FUNCTION SVF(SG5,SG6,SG7,SG8)
SVF=exp(SG5-1-log(SG5-SG6)-SG7*SG8)
END
此程序没有错误和warning了，但是就是运行不出结果，我是COMPAQ VISUAL FORTRAN.
回复此楼