Difference between revisions of "FORTRAN code for the Kolafa and Nezbeda equation of state"

From SklogWiki
Jump to: navigation, search
m (Removed white space.)
m (Added a comment re MAIN, and added spacers for column 6.)
Line 1: Line 1:
The following is the FORTRAN source code for the [[Lennard-Jones_equation_of_state#Kolafa and Nezbeda equation of state | Kolafa and Nezbeda equation of state]] for the [[Lennard-Jones model]]
+
The following are the FORTRAN source code functions (without a "main" program) for the [[Lennard-Jones_equation_of_state#Kolafa and Nezbeda equation of state | Kolafa and Nezbeda equation of state]] for the [[Lennard-Jones model]]
  
 
  c===================================================================
 
  c===================================================================
Line 18: Line 18:
 
       eta = PI/6.*rho * (dC(T))**3
 
       eta = PI/6.*rho * (dC(T))**3
 
       ALJ =  (dlog(rho)+betaAHS(eta)
 
       ALJ =  (dlog(rho)+betaAHS(eta)
    +  +rho*BC(T)/exp(gammaBH(T)*rho**2))*T
+
      +  +rho*BC(T)/exp(gammaBH(T)*rho**2))*T
    +  +DALJ(T,rho)
+
      +  +DALJ(T,rho)
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 28: Line 28:
 
       eta = PI/6. *rho*(dC(T))**3
 
       eta = PI/6. *rho*(dC(T))**3
 
       ALJres = (betaAHS(eta)
 
       ALJres = (betaAHS(eta)
    + +rho*BC(T)/exp(gammaBH(T)*Sqrt(rho)))*T
+
      + +rho*BC(T)/exp(gammaBH(T)*Sqrt(rho)))*T
    + +DALJ(T,rho)
+
      + +DALJ(T,rho)
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 38: Line 38:
 
       eta=PI/6. *rho*(dC(T))**3
 
       eta=PI/6. *rho*(dC(T))**3
 
       sum=((2.01546797*2+rho*(
 
       sum=((2.01546797*2+rho*(
    + (-28.17881636)*3+rho*(
+
      + (-28.17881636)*3+rho*(
    + 28.28313847*4+rho*
+
      + 28.28313847*4+rho*
    + (-10.42402873)*5)))
+
      + (-10.42402873)*5)))
    + +(-19.58371655*2+rho*(
+
      + +(-19.58371655*2+rho*(
    + +75.62340289*3+rho*(
+
      + +75.62340289*3+rho*(
    + (-120.70586598)*4+rho*(
+
      + (-120.70586598)*4+rho*(
    + +93.92740328*5+rho*
+
      + +93.92740328*5+rho*
    + (-27.37737354)*6))))/dsqrt(T)
+
      + (-27.37737354)*6))))/dsqrt(T)
    + + ((29.34470520*2+rho*(
+
      + + ((29.34470520*2+rho*(
    + (-112.35356937)*3+rho*(
+
      + (-112.35356937)*3+rho*(
    + +170.64908980*4+rho*(
+
      + +170.64908980*4+rho*(
    + (-123.06669187)*5+rho*
+
      + (-123.06669187)*5+rho*
    + 34.42288969*6))))+
+
      + 34.42288969*6))))+
    + ((-13.37031968)*2+rho*(
+
      + ((-13.37031968)*2+rho*(
    + 65.38059570*3+rho*(
+
      + 65.38059570*3+rho*(
    + (-115.09233113)*4+rho*(
+
      + (-115.09233113)*4+rho*(
    + 88.91973082*5+rho*
+
      + 88.91973082*5+rho*
    + (-25.62099890)*6))))/T)/T)*rho**2
+
      + (-25.62099890)*6))))/T)/T)*rho**2
 
         PLJ = ((zHS(eta)
 
         PLJ = ((zHS(eta)
    +  + BC(T)/exp(gammaBH(T)*rho**2)
+
      +  + BC(T)/exp(gammaBH(T)*rho**2)
    +  *rho*(1-2*gammaBH(T)*rho**2))*T
+
      +  *rho*(1-2*gammaBH(T)*rho**2))*T
    +  +sum )*rho
+
      +  +sum )*rho
 
         RETURN
 
         RETURN
 
         END
 
         END
Line 71: Line 71:
 
       eta=PI/6. *rho*d**3
 
       eta=PI/6. *rho*d**3
 
       sum= ((2.01546797+rho*(
 
       sum= ((2.01546797+rho*(
    + (-28.17881636)+rho*(
+
      + (-28.17881636)+rho*(
    + +28.28313847+rho*
+
      + +28.28313847+rho*
    + (-10.42402873))))
+
      + (-10.42402873))))
    + + (-19.58371655*1.5+rho*(
+
      + + (-19.58371655*1.5+rho*(
    + 75.62340289*1.5+rho*(
+
      + 75.62340289*1.5+rho*(
    + (-120.70586598)*1.5+rho*(
+
      + (-120.70586598)*1.5+rho*(
    + 93.92740328*1.5+rho*
+
      + 93.92740328*1.5+rho*
    + (-27.37737354)*1.5))))/dsqrt(T)
+
      + (-27.37737354)*1.5))))/dsqrt(T)
    + + ((29.34470520*2+rho*(
+
      + + ((29.34470520*2+rho*(
    + -112.35356937*2+rho*(
+
      + -112.35356937*2+rho*(
    +  170.64908980*2+rho*(
+
      +  170.64908980*2+rho*(
    + -123.06669187*2+rho*
+
      + -123.06669187*2+rho*
    + 34.42288969*2)))) +
+
      + 34.42288969*2)))) +
    + (-13.37031968*3+rho*(
+
      + (-13.37031968*3+rho*(
    +  65.38059570*3+rho*(
+
      +  65.38059570*3+rho*(
    +  -115.09233113*3+rho*(
+
      +  -115.09233113*3+rho*(
    + 88.91973082*3+rho*
+
      + 88.91973082*3+rho*
    + (-25.62099890)*3))))/T)/T) *rho*rho
+
      + (-25.62099890)*3))))/T)/T) *rho*rho
      ULJ = 3*(zHS(eta)-1)*dBHdT/d
+
      ULJ = 3*(zHS(eta)-1)*dBHdT/d
    + +rho*dB2BHdT/exp(gammaBH(T)*rho**2) +sum
+
      + +rho*dB2BHdT/exp(gammaBH(T)*rho**2) +sum
      RETURN
+
      RETURN
      END
+
      END
 
       DOUBLE PRECISION FUNCTION zHS(eta)
 
       DOUBLE PRECISION FUNCTION zHS(eta)
 
       implicit double precision (a-h,o-z)
 
       implicit double precision (a-h,o-z)
Line 101: Line 101:
 
       implicit double precision (a-h,o-z)
 
       implicit double precision (a-h,o-z)
 
       betaAHS = dlog(1-eta)/0.6
 
       betaAHS = dlog(1-eta)/0.6
    +  + eta*( (4.0/6*eta-33.0/6)*eta+34.0/6 ) /(1.-eta)**2
+
      +  + eta*( (4.0/6*eta-33.0/6)*eta+34.0/6 ) /(1.-eta)**2
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 110: Line 110:
 
       isT=1/dsqrt(T)
 
       isT=1/dsqrt(T)
 
       dLJ = ((( 0.011117524191338 *isT-0.076383859168060)
 
       dLJ = ((( 0.011117524191338 *isT-0.076383859168060)
    + *isT)*isT+0.000693129033539)/isT+1.080142247540047
+
      + *isT)*isT+0.000693129033539)/isT+1.080142247540047
    + +0.127841935018828*dlog(isT)
+
      + +0.127841935018828*dlog(isT)
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 118: Line 118:
 
       sT=dsqrt(T)
 
       sT=dsqrt(T)
 
       dC = -0.063920968*dlog(T)+0.011117524/T
 
       dC = -0.063920968*dlog(T)+0.011117524/T
    +    -0.076383859/sT+1.080142248+0.000693129*sT
+
      +    -0.076383859/sT+1.080142248+0.000693129*sT
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 125: Line 125:
 
       sT=dsqrt(T)
 
       sT=dsqrt(T)
 
       dCdT =  0.063920968*T+0.011117524+(-0.5*0.076383859
 
       dCdT =  0.063920968*T+0.011117524+(-0.5*0.076383859
    +  -0.5*0.000693129*T)*sT
+
      +  -0.5*0.000693129*T)*sT
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 133: Line 133:
 
       isT=1/dsqrt(T)
 
       isT=1/dsqrt(T)
 
       BC = (((((-0.58544978*isT+0.43102052)*isT
 
       BC = (((((-0.58544978*isT+0.43102052)*isT
    +  +.87361369)*isT-4.13749995)*isT+2.90616279)*isT
+
      +  +.87361369)*isT-4.13749995)*isT+2.90616279)*isT
    +  -7.02181962)/T+0.02459877
+
      +  -7.02181962)/T+0.02459877
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 142: Line 142:
 
       isT=1/dsqrt(T)
 
       isT=1/dsqrt(T)
 
       BCdT = ((((-0.58544978*3.5*isT+0.43102052*3)*isT
 
       BCdT = ((((-0.58544978*3.5*isT+0.43102052*3)*isT
    +  +0.87361369*2.5)*isT-4.13749995*2)*isT
+
      +  +0.87361369*2.5)*isT-4.13749995*2)*isT
    +  +2.90616279*1.5)*isT-7.02181962
+
      +  +2.90616279*1.5)*isT-7.02181962
 
       RETURN
 
       RETURN
 
       END
 
       END
Line 154: Line 154:
 
       implicit double precision (a-h,o-z)
 
       implicit double precision (a-h,o-z)
 
       DALJ = ((+2.01546797+rho*(-28.17881636
 
       DALJ = ((+2.01546797+rho*(-28.17881636
    + +rho*(+28.28313847+rho*(-10.42402873))))
+
      + +rho*(+28.28313847+rho*(-10.42402873))))
    + +(-19.58371655+rho*(75.62340289+rho*((-120.70586598)
+
      + +(-19.58371655+rho*(75.62340289+rho*((-120.70586598)
    + +rho*(93.92740328+rho*(-27.37737354)))))/dsqrt(T)
+
      + +rho*(93.92740328+rho*(-27.37737354)))))/dsqrt(T)
    + + ( (29.34470520+rho*((-112.35356937)
+
      + + ( (29.34470520+rho*((-112.35356937)
    + +rho*(+170.64908980+rho*((-123.06669187)
+
      + +rho*(+170.64908980+rho*((-123.06669187)
    + +rho*34.42288969))))
+
      + +rho*34.42288969))))
    + +(-13.37031968+rho*(65.38059570+
+
      + +(-13.37031968+rho*(65.38059570+
    + rho*((-115.09233113)+rho*(88.91973082
+
      + rho*((-115.09233113)+rho*(88.91973082
    + +rho* (-25.62099890)))))/T)/T) *rho*rho
+
      + +rho* (-25.62099890)))))/T)/T) *rho*rho
 
       RETURN
 
       RETURN
 
       END
 
       END

Revision as of 12:13, 25 November 2010

The following are the FORTRAN source code functions (without a "main" program) for the Kolafa and Nezbeda equation of state for the Lennard-Jones model

c===================================================================
C      Package supplying the thermodynamic properties of the
C      LENNARD-JONES fluid
c
c      J. Kolafa, I. Nezbeda, Fluid Phase Equil. 100 (1994), 1
c
c      ALJ(T,rho)...Helmholtz free energy (including the ideal term)
c      PLJ(T,rho)...Pressure
c      ULJ(T,rho)...Internal energy
c===================================================================
      DOUBLE PRECISION FUNCTION ALJ(T,rho)
C      Helmholtz free energy (including the ideal term)
c
      implicit double precision (a-h,o-z)
      data pi /3.141592654d0/
      eta = PI/6.*rho * (dC(T))**3
      ALJ =  (dlog(rho)+betaAHS(eta)
     +  +rho*BC(T)/exp(gammaBH(T)*rho**2))*T
     +  +DALJ(T,rho)
      RETURN
      END
C/* Helmholtz free energy (without ideal term) */
      DOUBLE PRECISION FUNCTION ALJres(T,rho)
      implicit double precision (a-h,o-z)
      data pi /3.141592654d0/
      eta = PI/6. *rho*(dC(T))**3
      ALJres = (betaAHS(eta)
     + +rho*BC(T)/exp(gammaBH(T)*Sqrt(rho)))*T
     + +DALJ(T,rho)
      RETURN
      END
C/* pressure */
      DOUBLE PRECISION FUNCTION PLJ(T,rho)
      implicit double precision (a-h,o-z)
      data pi /3.141592654d0/
      eta=PI/6. *rho*(dC(T))**3
      sum=((2.01546797*2+rho*(
     + (-28.17881636)*3+rho*(
     + 28.28313847*4+rho*
     + (-10.42402873)*5)))
     + +(-19.58371655*2+rho*(
     + +75.62340289*3+rho*(
     + (-120.70586598)*4+rho*(
     + +93.92740328*5+rho*
     + (-27.37737354)*6))))/dsqrt(T)
     + + ((29.34470520*2+rho*(
     + (-112.35356937)*3+rho*(
     + +170.64908980*4+rho*(
     + (-123.06669187)*5+rho*
     + 34.42288969*6))))+
     + ((-13.37031968)*2+rho*(
     + 65.38059570*3+rho*(
     + (-115.09233113)*4+rho*(
     + 88.91973082*5+rho*
     + (-25.62099890)*6))))/T)/T)*rho**2
       PLJ = ((zHS(eta)
     +  + BC(T)/exp(gammaBH(T)*rho**2)
     +  *rho*(1-2*gammaBH(T)*rho**2))*T
     +  +sum )*rho
       RETURN
       END
C/* internal energy */
      DOUBLE PRECISION FUNCTION ULJ( T, rho)
      implicit double precision (a-h,o-z)
       data pi /3.141592654d0/
      dBHdT=dCdT(T)
      dB2BHdT=BCdT(T)
      d=dC(T)
      eta=PI/6. *rho*d**3
      sum= ((2.01546797+rho*(
     + (-28.17881636)+rho*(
     + +28.28313847+rho*
     + (-10.42402873))))
     + + (-19.58371655*1.5+rho*(
     + 75.62340289*1.5+rho*(
     + (-120.70586598)*1.5+rho*(
     + 93.92740328*1.5+rho*
     + (-27.37737354)*1.5))))/dsqrt(T)
     + + ((29.34470520*2+rho*(
     + -112.35356937*2+rho*(
     +  170.64908980*2+rho*(
     + -123.06669187*2+rho*
     + 34.42288969*2)))) +
     + (-13.37031968*3+rho*(
     +  65.38059570*3+rho*(
     +  -115.09233113*3+rho*(
     + 88.91973082*3+rho*
     + (-25.62099890)*3))))/T)/T) *rho*rho
      ULJ = 3*(zHS(eta)-1)*dBHdT/d
     + +rho*dB2BHdT/exp(gammaBH(T)*rho**2) +sum
      RETURN
      END
      DOUBLE PRECISION FUNCTION zHS(eta)
      implicit double precision (a-h,o-z)
      zHS = (1+eta*(1+eta*(1-eta/1.5*(1+eta)))) / (1-eta)**3
      RETURN
      END
      DOUBLE PRECISION FUNCTION betaAHS( eta )
      implicit double precision (a-h,o-z)
      betaAHS = dlog(1-eta)/0.6
     +  + eta*( (4.0/6*eta-33.0/6)*eta+34.0/6 ) /(1.-eta)**2
      RETURN
      END
C /* hBH diameter */
      DOUBLE PRECISION FUNCTION dLJ(T)
      implicit double precision (a-h,o-z)
      DOUBLE PRECISION IST
      isT=1/dsqrt(T)
      dLJ = ((( 0.011117524191338 *isT-0.076383859168060)
     + *isT)*isT+0.000693129033539)/isT+1.080142247540047
     + +0.127841935018828*dlog(isT)
      RETURN
      END
      DOUBLE PRECISION FUNCTION dC(T)
      implicit double precision (a-h,o-z)
      sT=dsqrt(T)
      dC = -0.063920968*dlog(T)+0.011117524/T
     +     -0.076383859/sT+1.080142248+0.000693129*sT
      RETURN
      END
      DOUBLE PRECISION FUNCTION dCdT( T)
      implicit double precision (a-h,o-z)
      sT=dsqrt(T)
      dCdT =   0.063920968*T+0.011117524+(-0.5*0.076383859
     +   -0.5*0.000693129*T)*sT
      RETURN
      END
      DOUBLE PRECISION FUNCTION BC( T)
      implicit double precision (a-h,o-z)
      DOUBLE PRECISION isT
      isT=1/dsqrt(T)
      BC = (((((-0.58544978*isT+0.43102052)*isT
     +  +.87361369)*isT-4.13749995)*isT+2.90616279)*isT
     +  -7.02181962)/T+0.02459877
      RETURN
      END
      DOUBLE PRECISION FUNCTION BCdT( T)
      implicit double precision (a-h,o-z)
      DOUBLE PRECISION iST
      isT=1/dsqrt(T)
      BCdT = ((((-0.58544978*3.5*isT+0.43102052*3)*isT
     +  +0.87361369*2.5)*isT-4.13749995*2)*isT
     +  +2.90616279*1.5)*isT-7.02181962
      RETURN
      END
      DOUBLE PRECISION FUNCTION gammaBH(X)
      implicit double precision (a-h,o-z)
      gammaBH=1.92907278
      RETURN
      END
      DOUBLE PRECISION FUNCTION DALJ(T,rho)
      implicit double precision (a-h,o-z)
      DALJ = ((+2.01546797+rho*(-28.17881636
     + +rho*(+28.28313847+rho*(-10.42402873))))
     + +(-19.58371655+rho*(75.62340289+rho*((-120.70586598)
     + +rho*(93.92740328+rho*(-27.37737354)))))/dsqrt(T)
     + + ( (29.34470520+rho*((-112.35356937)
     + +rho*(+170.64908980+rho*((-123.06669187)
     + +rho*34.42288969))))
     + +(-13.37031968+rho*(65.38059570+
     + rho*((-115.09233113)+rho*(88.91973082
     + +rho* (-25.62099890)))))/T)/T) *rho*rho
      RETURN
      END


Reference

Source.png This page contains computer source code. If you intend to compile and use this code you must check for yourself the validity of the code. Please read the SklogWiki disclaimer.
40px-Stop hand nuvola.svg.png This page contains numerical values and/or equations. If you intend to use ANY of the numbers or equations found in SklogWiki in any way, you MUST take them from the original published article or book, and cite the relevant source accordingly.