[Fortran] 纯文本查看 复制代码
program enthalpy
implicit none
!This program aim to calculate the enthalpy of h2s hydrate dissociation in the presence of water
double precision P,T,xh2s,n
double precision dPdT,H1,V,Hdiss,Hsolve
print*,"please imput temperature (K)"
read*,T
print*,"please imput pressure (MPa)"
read*,P
print*,"please imput hydation number"
read*,n
print*,"please imput mole fraction of h2s in liquid"
read*,xh2s
call dissenthalpy(Hdiss,T,dPdT,V)
call solveenthalpy(Hsolve,xh2s,n)
H1=Hdiss+Hsolve
print*,"H1=",H1
pause
! return
end
subroutine Derivative(dPDT,T)
double precision dPdT,T,e(4)
data e/-436.1,4.6,-0.018,0.0000196/
dPdT=e(1)+e(2)*T+e(3)*T**2+e(4)*T**3 !MPa/K
print*,"dPdT=",dPdT
return
end
subroutine Totalmolarvolume(V,Vliquid,Vgas,Vhydrate)
double precision V,Vliquid,Vgas,Vhydrate,xh2s,n,parameter1,parameter2,density2,l(4) !density1,molecularweight2,V4,molecularweight,h(6)
!data h/44.202,2544.564,-35.957,0.191,-4.488E-4,3.965E-7/
data l/21.90583,-0.21947,7.33695E-4,-8.18347E-7/
!calculating the molar volume change of solution
parameter1=0.0000045
parameter2=0.000018
Vliquid=xh2s*parameter1+parameter2 !m3/mol
!calculating the molar volume change of hydrate
! density1=1050 !kg/m3 !this value is guess
! molecularweight2=0.1461 !kg/mol
Vhydrate=1.391E-4 !m3/mol, Vhydrate=molecularweight/density1
!calculating the molar volume change of vapor
!density2=h(1)+h(2)*T+h(3)*T**2+h(4)*T**3+h(5)*T**4+h(6)*T**5 !kg/m3
!molecularweight=0.034 !kg/mol
! Vgas=molecularweight/density2 !m3/mol
Vgas=l(1)+l(2)*T+l(3)*T**2+l(4)*T**3
!calculating the total molar volume
V=(1-(n*xh2s)/(1-xh2s))*Vgas+n*Vliquid-Vhydrate !m3/mol
print*,"Vliquid=",Vliquid
print*,"Vhydrate=",Vhydrate
print*,"Vgas=",Vgas
print*,"V=",V
return
end
subroutine dissenthalpy(Hdiss,T,dPdT,V)
double precision Hdiss,T,dPdT,V,Vliquid,Vgas,Vhydrate
call Derivative(dPDT,T)
call Totalmolarvolume(V,Vliquid,Vgas,Vhydrate)
Hdiss=1000000.0*dPdT*T*V
print*,"Hdiss=",Hdiss
return
end
subroutine solveenthalpy(Hsolve,xh2s,n)
double precision Hsolve,xh2s,n
Hsolve=17960.0*n*xh2s/(1-xh2s) !J/mol
print*,"Hsolve=",Hsolve
return
end