Question

Butane liquid and vapor coexist at 370.0K and 14.35 bar. The densities of the liquid and vapor phases are 8.128 mol-L^-1 and 0.6313 mol-L^-1, respectively. Use the van der Waals equation, the Redlich-Kwong equation, and the Peng-Robinson equation to calculate these densities. Take alpha = 16.44 bar-L^2-mol^-2 and beta = 0.07245 L-mol^-1 for the Peng-Robinson equation.

(Chapter 16, problem 24 in the Physical Chemistry, a molecular approach book)

Answer 1

Program in scientific python:

# -*- coding: utf-8 -*-

"""

Created on Mon Oct 05 18:01:59 2015

@author: Sebastian

a) Van der Waals

b) Peng Robinson

c) Soave RK

"""

from scipy.constants import R as Rgi

from numpy import *

from scipy.optimize import fsolve

"""systems conditions"""

T = 370 #k

P = 1.435E6 #Pa

"""Critical conditions of butane"""#Cheric data base

Tc = 4.25120E+02 #K

Pc = 3.796000E+006 #Pa

omega = 1.99000E-01

def FunZVdWvap(T,P):

global r

a = (27*Rgi**2*Tc**2)/(64*Pc)

b = (Rgi*Tc)/(8*Pc)

A = (P*a)/((Rgi*T)**2)

B = (P*b)/(Rgi*T)

beta2 = -(B+1)

beta1 = A

beta0 = -A*B

pol = [1,beta2,beta1,beta0]

r = roots(pol)

Z = real(r[isreal(r)])[0]

return Z

def FunZVdWliq(T,P):

global r

a = (27*Rgi**2*Tc**2)/(64*Pc)

b = (Rgi*Tc)/(8*Pc)

A = (P*a)/((Rgi*T)**2)

B = (P*b)/(Rgi*T)

beta2 = -(B+1)

beta1 = A

beta0 = -A*B

pol = [1,beta2,beta1,beta0]

r = roots(pol)

Z = real(r[isreal(r)])[2]

return Z

Denvap = P/(Rgi*T*FunZVdWvap(T,P))*(1.0/1000)

Denliq = P/(Rgi*T*FunZVdWliq(T,P))*(1.0/1000)

print '{0}'.format('Van der Waals Equation')

print '{0}'.format('Vapor density (mol/L), Liquid density(mol/L) ')

print '{0} {1}'.format(Denvap,Denliq)

"""Peng Robinson EoS"""

def PR_EOSvap(T,P):

b = 0.07780*(Rgi*Tc)/(Pc)

kappa = 0.37464 + 1.54226*omega -0.26992*omega**2

alpha = (1 + kappa*(1-sqrt(T/Tc)))**2

a = 0.45724*alpha*((Rgi*Tc)**2)/Pc

B = (b*P)/(Rgi*T)

A = (a*P)/(Rgi*T)**2

beta2 = B-1

beta1 = A - 3*(B**2) -2*B

beta0 = B**3 + B**2 - A*B

pol = [1,beta2,beta1,beta0]

r = roots(pol)

Z = real(r[isreal(r)])[0] #plot for vapor

return Z

def PR_EOSliq(T,P):

b = 0.07780*(Rgi*Tc)/(Pc)

kappa = 0.37464 + 1.54226*omega -0.26992*omega**2

alpha = (1 + kappa*(1-sqrt(T/Tc)))**2

a = 0.45724*alpha*((Rgi*Tc)**2)/Pc

B = (b*P)/(Rgi*T)

A = (a*P)/(Rgi*T)**2

beta2 = B-1

beta1 = A - 3*(B**2) -2*B

beta0 = B**3 + B**2 - A*B

pol = [1,beta2,beta1,beta0]

r = roots(pol)

Z = real(r[isreal(r)])[2] #plot for liquid

return Z

Denvap2 = P/(Rgi*T*PR_EOSvap(T,P))*(1.0/1000)

Denliq2 = P/(Rgi*T*PR_EOSliq(T,P))*(1.0/1000)

print '{0}'.format('Peng Robinson EoS')

print '{0}'.format('Vapor density (mol/L), Liquid density(mol/L) ')

print '{0} {1}'.format(Denvap2,Denliq2)

"""SRK EoS"""

def SRK_EOSvap(T,P):

b = 0.07780*(Rgi*Tc)/(Pc)

kappa = 0.37464 + 1.54226*omega -0.26992*omega**2

alpha = (1 + kappa*(1-sqrt(T/Tc)))**2

a = 0.45724*alpha*((Rgi*Tc)**2)/Pc

B = (b*P)/(Rgi*T)

A = (a*P)/(Rgi*T)**2

beta2 = -1

beta1 = A - B - B**2

beta0 = -A*B

pol = [1,beta2,beta1,beta0]

r = roots(pol)

Z = real(r[isreal(r)])[0] #plot for vapor

return Z

def SRK_EOSliq(T,P):

b = 0.07780*(Rgi*Tc)/(Pc)

kappa = 0.37464 + 1.54226*omega -0.26992*omega**2

alpha = (1 + kappa*(1-sqrt(T/Tc)))**2

a = 0.45724*alpha*((Rgi*Tc)**2)/Pc

B = (b*P)/(Rgi*T)

A = (a*P)/(Rgi*T)**2

beta2 = -1

beta1 = A - B - B**2

beta0 = -A*B

pol = [1,beta2,beta1,beta0]

r = roots(pol)

Z = real(r[isreal(r)])[2] #plot for liquid

return Z

Denvap3 = P/(Rgi*T*SRK_EOSvap(T,P))*(1.0/1000)

Denliq3 = P/(Rgi*T*SRK_EOSliq(T,P))*(1.0/1000)

print '{0}'.format('Soave RK EoS')

print '{0}'.format('Vapor density (mol/L), Liquid density(mol/L) ')

print '{0} {1}'.format(Denvap3,Denliq3)

Program outputs:

Van der Waals Equation
Vapor density (mol/L), Liquid density(mol/L)
0.574120847444 4.78587403725
Peng Robinson EoS
Vapor density (mol/L), Liquid density(mol/L)
0.631976283254 8.11360467905
Soave RK EoS
Vapor density (mol/L), Liquid density(mol/L)
0.649763368867 9.20326108369
>>>