问题：srk ，pr 方程闪蒸计算 fortran源代码
类型：交流
提问：shx8282
等级：▲
版块：化工工艺设计(wanghong,)
信誉：100%
回复：2
阅读：905
时间：2009-11-09 07:49:47 编辑加入/取消收藏订制/取消短消息举报该贴

! flash2.f90
!
! FUNCTIONS:
! flash2 - Entry point of console application.
!
! Example of displaying 'Hello World' at execution time.
!

!****************************************************************************
!
! PROGRAM: flash2
!
! PURPOSE: Entry point for 'Hello World' sample console application.
!
!****************************************************************************

!srk ，pr 方程闪蒸计算 fortran源代码

program flash2
chArActer*100 zhushi
chArActer*3 xh
REAL,AllocAtAble :: PC(:),TC(:),WM(:),OMEGA(:),X(:),Y(:),Z(:),ak(:,:),PHIL(:) ,PHIV(:),Bk(:),Y1(:),sk(:)
open (22,file='sr.txt')
READ(22,*) zhushi
READ(22,*)n
READ(22,*) zhushi
READ(22,*)p
READ(22,*) zhushi
READ(22,*)t
AllocAte( PC(n),TC(n),OMEGA(n),X(n),Y(n),Z(n),ak(n,n),PHIL(n) ,PHIV(n),Bk(n),Y1(n),sk(n),WM(n))
READ(22,*) zhushi
do i=1 ,n
READ(22,*)pc(i),tc(i),OMEGA(i),Z(i),WM(I)
enddo
READ(22,*) zhushi
READ(22,*)ak(:,:)
READ(22,*) zhushi
READ(22,*)neos
close(22)

call paodianpaodian(n,p,PC,TC,OMEGA,z,AK)!(n,p,t,PC,TC,OMEGA,z,AK,x,y,e,PHIL ,PHIV,sk,yam,ybm,yaa,ybb,qam,qbm,qaa,qbb,neos,yz,qz)
!call flash(n,p,t,PC,TC,OMEGA,z,AK,x,y,e)
open (22,file='jg.txt')
if (e==-1.0)then
write(22,*)'闪蒸计算失败，处于单相区'
stop
end if
if(neos==0) then
write(22,*) 'PR方程'
else
write(22,*) 'SRK方程'
ENDIF
do i=1,n
zppc=qppc+z(I)*pc(i)
zttc=qttc+z(I)*tc(i)
zomg=omg+z(I)*OMEGA(i)
zWM=zWM+z(I)*WM(I)
enddo

write(22,*) ' 温度(K ) ',t
write(22,*) ' 压力(PA) ',p
write(22,*) ' 气化率 ',e
write(22,*) ' 临界压力(PA) ',zppc
write(22,*) ' 临界温度(K ) ',zttc
write(22,*) ' 偏心因子 ',zomg
write(22,*) ' 分子量 ',zWM

write(22,*)'-------------------------------------------------------------------------'
write(22,*)'-------------------------------------------------------------------------'
write(22,*)' '
yppc=0.
qppc=0.
yttc=0.
qttc=0.
qomg=0.
yomg=0.
do i=1,n
yppc=yppc+x(i)*pc(i)
qppc=qppc+y(i)*pc(i)
yttc=yttc+x(i)*tc(i)
qttc=qttc+y(i)*tc(i)
qomg=omg+y(i)*OMEGA(i)
yomg=omg+x(i)*OMEGA(i)
YWM=YWM+X(I)*WM(I)
QWM=QWM+Y(I)*WM(I)
enddo
ymv=yz *8.314*t /p
qmv=qz *8.314*t /p
ymidu=YWM/ymv/1000
qmidu=qWM/qmv/1000
write(22,*) ' 液相 ' ,'气相 '
write(22,*)'-------------------------------------------------------------------------'

write(22,*) ' 摩尔分率 ',1-e, ' ',e
write(22,*) ' 临界压力(PA) ' ,yppc,' ',qppc
write(22,*) ' 临界温度(K ) ' ,yttc,' ',qttc
write(22,*) ' 偏心因子 ' ,yomg,' ',qomg
write(22,*) ' 压缩因子 ' ,yz,' ',qz
write(22,*) ' 分子量 ' ,YWM,' ',QWM
write(22,*)' 密度(kg/m3)' ,ymidu,' ',qmidu
write(22,*) ' 摩尔体积(m3) ' ,ymv,' ',qmv
write(22,*) ' 方程am ',yam, ' ',qam
write(22,*) ' 方程bm ',ybm, ' ',qbm
write(22,*) ' 方程aa ',yaa, ' ',qaa
write(22,*) ' 方程bb ',ybb, ' ',qbb
write(22,*)''
write(22,*) ' 组成液相 ' ,'气相 '
write(22,*)'-------------------------------------------------------------------------'
do i=1,n
write(xh,'(i3)')i
write(22,*) '组分'//xh ,' ',x(i),' ',y(i)!,sk(i),PHIV(i),PHIL(i)
end do
write(22,*)''
write(22,*) ' 逸度液相 ' ,'气相 '
write(22,*)'-------------------------------------------------------------------------'
do i=1,n
write(22,*) ' ',PHIL(i),' ',PHIV(i)!,sk(i),PHIV(i),PHIL(i)
end do
write(22,*)''
write(22,*) ' 平衡常数液相 ' ,'气相 '
write(22,*)'-------------------------------------------------------------------------'

do i=1,n
write(22,*)' ',sk(i),' ',sk(i)!,sk(i),PHIV(i),PHIL(i)
end do
close(22)
end

SUBROUTINE flash(n,p,t,PC,TC,OMEGA,z,AK,x,y,e,PHIL ,PHIV,sk,yam,ybm,yaa,ybb,qam,qbm,qaa,qbb,neos,yz,qz)
real:: PC(n),TC(n),OMEGA(n),X(n),Y(n),Z(n),AK(n,n),PHIL(n) ,PHIV(n),Bk(n),Y1(n),sk(n)
open (23,file='gc.txt')
20 do i = 1 ,n
tr=t/tc(i)
pr=p/pc(i)
sk(i) = Exp(5.37 * (1 + OMEGA(i)) * (1 - 1/tr)) /pr
end do
nkz=0
15 continue
e=0.6
9 kjs=0
10 continue
kjs=kjs+1
f=0
df=0
sum=0
sum1=0
do i=1,n
smm=(sk(i)-1)/(1+(sk(i)-1)*e)
f=f+z(i)*smm
df=df-z(i)*smm**2
sum=sum +sk(i)*z(i)
sum1=sum1 +z(i)/sk(i)
write(23,*)sk(i)
end do
write(23,*)'------------------------------'
write(23,*)e
write(23,*)'------------------------------'
e1=e
if(df==0. .or. kjs>100) goto 100
e=e-f/df
if (abs(e-e1)/e>0.0001)goto 10
write(23,*)'=========================='
sum1=0
sum2=0
if(e==e2.and. e==e3)goto 30
e2=e
e3=e2
do i=1,n
x(i)=z(i)/(1+(sk(i)-1)*e)
y(i)=sk(i)*x(i)
sum1=sum1+x(i)
sum2=sum2+y(i)
end do
do i=1,n
x(i)=x(i)/sum1
y(i)=y(i)/sum2
end do
zz=0.9
if (neos==1)then
CALL srkEOS(n,T,P,TC,PC,OMEGA,AK,y,zz,PHIV,qam,qbm,qaa,qbb,qz)
else
CALL PREOS(n,T,P,TC,PC,OMEGA,AK,y,zz,PHIV,qam,qbm,qaa,qbb,qz)
end if
z1=zz
zz=0.1
if (neos==1)then
CALL srkEOS(n,T,P,TC,PC,OMEGA,AK,X,zz,PHIL,yam,ybm,yaa,ybb,yz)
else
CALL PREOS(n,T,P,TC,PC,OMEGA,AK,x,zz,PHIl,yam,ybm,yaa,ybb,yz)
end if
z2=zz
do i=1,n
sK(I)=PHIL(I)*y(i)/PHIV(I)/x(i)
end do
mtc=0
do i=1,n
if(PHIL(I)-PHIV(I)>0.01) mtc=1
end do
if (mtc==0)goto 30
nkz=nkz+1
if(nkz>200)then
e=-1
return
end if
goto 15
30 continue
return
100 e=0
end
SUBROUTINE PREOS(N,T,P,TC,PC,OMEGA, AK,Y,Z,PHI,am,bm,aa,bb,zz)
real:: TC(n),PC(n),OMEGA(n),AK(n,n),PHI(n),Y(n),A(n),B(n)
real zzz(3)
R=8.314
! Calculate Equation Parameter a and b
sum=0
DO 10 I=1, N
TR=T/TC(I)
AM=0.37464+1.54226 * OMEGA(I)-0.26992 * OMEGA(I) ** 2
! AM1= 0.3796 + 1.485 * OMEGA(i)- 0.1644 * OMEGA(i)**2+ 0.01667 * OMEGA(i)**3
ALPHA = (1.0+AM * (1.0-SQRT(TR))) ** 2
A(I) = ALPHA * 0.45724 * R* R *TC(I) * TC(I)/PC(I)
B(I) = 0.07780*R*TC(I)/PC(I)
sum=sum +tc(i)*y(i)
10 CONTINUE
AM = 0.0
BM=0.0
DO 20 I = 1,N
BM = BM+Y(I)*B(I)
DO 20 J=1,N
AM =AM+ Y(I) * Y(J) * (1.0-AK(I,J)) * SQRT(A(I) * A(J))
20 CONTINUE
! Calculate Compressibility Factor 2

AA=AM * P/(R *T) ** 2
BB=BM*P/(R*T)
cAll ROOT1(-(1. -BB) ,(AA -2. * BB - 3. * BB * BB) &
, -(AA * BB-BB ** 2-BB**3),xz,Zzz,DT)
if (xz==0.)then
if (z>0.5)then
xz =zzz(3)
else
xz=zzz(1)
end if
end if
zz=xz
V =ZZ*R*T/P
z=zz
! Calculate Fugacity Coefficients PHI CD
DO 40 I=1,N
BBI=B(I)*P/(R*T)
PH1 = BBI/BB * (ZZ- 1.0)- ALOG(ZZ - BB)
SUMYA = 0.0
DO 50 J = 1,N
SUMYA=SUMYA +Y(J) * (1-AK(I,J)) * SQRT(A(I) * A(J))
50 CONTINUE
BSUM=BBI/BB-2. *SUMYA/AM
Q=SQRT(2.0)
PHI(I)=PH1+AA/(BB * 2. * Q) * ALOG( (ZZ+ (1. + Q) * BB)/(Zz+ (1. - Q) * BB))*BSUM
PHI(I)=EXP(PHI(I))*y(i)*p
40 CONTINUE
RETURN
END
! ****************************************************************************************** !
! *************** 程序: ! pc !**************** !
! *************** 参数: 1=4540000*1.01325 !**************** !
! ****************************************************************************************** !
SUBROUTINE srkEOS(N,T,P,TC,PC,OMEGA, AK,Y,Z,PHI,am,bm,aa,bb,zz)
real:: TC(n),PC(n),OMEGA(n),AK(n,n),PHI(n),Y(n),A(n),B(n)
real zzz(3)
R=8314
! Calculate Equation Parameter a and b
sum=0
DO 10 I=1, N
TR=T/TC(I)
AM=0.48508+1.55171 * OMEGA(I)-0.15613 * OMEGA(I) ** 2
! AM1= 0.3796 + 1.485 * OMEGA(i)- 0.1644 * OMEGA(i)**2+ 0.01667 * OMEGA(i)**3
ALPHA = (1.0+AM * (1.0-SQRT(TR))) ** 2
A(I) = ALPHA * 0.42728 * R* R *TC(I) * TC(I)/PC(I)
B(I) = 0.08664*R*TC(I)/PC(I)
sum=sum +tc(i)*y(i)
10 CONTINUE
AM = 0.0
BM=0.0
DO 20 I = 1,N
BM = BM+Y(I)*B(I)
DO 20 J=1,N
AM =AM+ Y(I) * Y(J) * (1.0-AK(I,J)) * SQRT(A(I) * A(J))
20 CONTINUE
! Calculate Compressibility Factor 2

AA=AM * P/(R *T) ** 2
BB=BM*P/(R*T)
cAll ROOT1(-1.0 ,(AA - BB - BB * BB) &
, -(AA * BB),xz,Zzz,DT)
if (xz==0.)then
if (z>0.5)then
xz =zzz(3)
else
xz=zzz(1)
end if
end if
zz=xz
V =ZZ*R*T/P
z=zz
! Calculate Fugacity Coefficients PHI CD
DO 40 I=1,N
BBI=B(I)
PH1 = BBI/Bm * (ZZ- 1.0)- ALOG(ZZ - BB)
SUMYA = 0.0
DO 50 J = 1,N
SUMYA=SUMYA +Y(J) * (1-AK(I,J)) * SQRT(A(I) * A(J))
50 CONTINUE
BSUM=BBI/Bm-2. *SUMYA/AM
Q=SQRT(2.0)
PHI(I)=PH1+AA/(BB ) * ALOG( 1+bb/zz)*BSUM
PHI(I)=EXP(PHI(I))*y(i)*p
40 CONTINUE
RETURN
END

SUBROUTINE ROOT1(P,Q,R,ZZ,Z,DT)
REAL M,N,N1,N2,Z(3)
DATA PI /3.14159/
CU(X)=X/ABS(X)*ABS(X)**0.333333
P2=PI/2.0
P3=PI/3.0
M=Q-P*P/3.
N=R-P*Q/3.+2.*P**3/27.
DT=N*N/4.-(-M)**3/27.
IF (DT) 100,100,110
100 FX=-N/SQRT((-M)**3/27.)/2.
FX=-ATAN(FX/SQRT(-FX*FX+1.))+P2
FX=FX/3.
zM=2.*SQRT(-M/3.)
DO 101 i =1,3
ZZ=ZM*COS(FX+2.*I*P3)
101 z(I)=ZZ-P/3.
DO 103 i=1,2
DO 102 J=I+1,3
If ( z(i) .GT. Z(J) ) THEN
ZZ=Z(J)
Z(J)=Z(I)
z(I)=Zz
END IF
102 CONTINUE
103 CONTINUE
zz=0.
RETURN
110 ZZ=SQRT(DT)
N1=-N/2.+ZZ
n2=-N/2.-ZZ
IF ( ABS(N2) .LT. 1.E-12 ) N2=1.E-12
ZZ= CU(N1) + CU(N2)
ZZ= ZZ-P/3.
DO 106 i=1,3
106 Z(I)=0.0
RETURN
END
来自 http://www.shihua100.com

[该帖子已被shx8282在2010-01-03 10:37:18编辑过]

附件：烃类闪蒸计算工具.rar (下载一次扣除可用分1分。附件提供人得到0.5分)

回复人：452612957,▲ (应化跨化工的考研学子) 时间：2009-11-27 11:35:39 编辑	1楼

回复人：452612957,▲ (应化跨化工的考研学子) 时间：2009-11-27 11:36:53 编辑	2楼
很好，不错谢谢楼主！

问题讨论没有结束...

您尚未进入本论坛,登录之后才可以回贴

7msec