from gekko import GEKKO
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.optimize import fsolve
from scipy.optimize import curve_fit
from scipy.interpolate import interp1d
#from numpy import *
R = 1.9858775 # cal/mol K
def vaporP(T: float, component: list):
'''
:param T: temp [K]
:param component: [vapor pressure parameters]
:return: Vapor Pressure [Pa]
'''
A = component[0]
B = component[1]
C = component[2]
D = component[3]
E = component[4]
return np.exp(A + B/T + C*np.log(T) + D*T**E)
def liquidMolarVolume(T: float, component: list):
'''
:param T: temp in [K}
:param component: component liquid molar volume
: constants in list or array
:return: liquid molar volume []
'''
A = component[0]
B = component[1]
C = component[2]
D = component[3]
return 1 / A * B ** (1 + (1 - T / C) ** D)
def Lij(aij: float, T: float, Vlvali: float,\
Vlvalj: float):
# Vl1 = liquidMolarVolume(T, Vlvali)
# Vl2 = liquidMolarVolume(T, Vlvalj)
Vl1 = Vlvali
Vl2 = Vlvalj
return Vl2 / Vl1 * np.exp(-aij / (R * T))
def Gamma(a12: float, a21: float, x1: float,\
Temp: float, Component: int):
'''
:param Component:
:param a12:
:param a21:
:param x1:
:param Temp:
:return:
'''
# Liquid molar volume
Vethanol = 5.8492e-2
VCyclochexane = 0.10882
x2 = 1.0 - x1
L12 = Lij(a12, Temp, Vethanol, VCyclochexane)
L21 = Lij(a21, Temp, VCyclochexane, Vethanol)
L11 = 1.0
L22 = 1.0
if Component == 1:
A = x1 + x2 * L12
B = x1 * L11 / (x1 * L11 + x2 * L12) \
+x2 * L21 / (x1 * L21 + x2 * L22)
elif Component == 2:
A = x2 + x1 * L21
B = x2 * L22 / (x2 * L22 + x1 * L21) \
+x1 * L12 / (x2 * L12 + x1 * L11)
return np.exp(1.0 - np.log(A) - B)
# yi * P = xi * γi * Pisat
# Component 1: Ethanol
# Component 2: Cyclohexane
def massTOMoleF(x: list):
#Converts RI to mole fraction of Ethanol
MW_Ethanol = 46.06844
MW_cyclohexane = 84.15948
xeth = 49.261*x**2 - 152.51*x + 117.32
xcyc = 1-xeth
moleth = xeth / MW_Ethanol
molcyc = xcyc/ MW_cyclohexane
return moleth/ (moleth + molcyc)
#Read in actual data
#filename = 'RealData.csv'
filename = 'http://apmonitor.com/me575/uploads/Main/vle_data.txt'
data = pd.read_csv(filename)
temp_data = data['Temp_D'].values + 273.15 # K
IndexMeasuredY = data['Distills_RI']
IndexMeasuredX = data['Bottoms_RI']
x1_data = massTOMoleF(IndexMeasuredX)
x1_data[0] = 0
x1_data[1] = 1
y1_data = massTOMoleF(IndexMeasuredY)
y1_data[0] = 0
y1_data[1] = 1
# Ethanol & Cyclohexane Properties ---------
# Psat paramater values
PsatC2H5OH = [73.304, -7122.3, -7.1424, \
2.8853e-6, 2.0]
PsatC6H12 = [51.087, -5226.4, -4.2278, \
9.7554e-18, 6.0]
# thermo Liquid Molar Volume paramaters
Vl_C2H5OH = [1.6288, 0.27469, 514.0, 0.23178]
Vl_C6H12 = [0.88998, 0.27376, 553.8, 0.28571]
# Build Model ------------------------------
m = GEKKO()
# Define Variables & Parameters
a12_ = m.FV(value = 2026.3)
a12_.STATUS = 1
a21_ = m.FV(value = 390.4)
a21_.STATUS = 1
Temp = m.Var(333)
x1 = m.Param(value = np.array(x1_data))
y1 = m.CV(value = np.array(y1_data))
y1.FSTATUS = 1 # min fstatus * (meas-pred)^2
y2 = m.Intermediate(1.0 - y1)
P = 85900 # Pa in Provo, UT
# Calc Vapor pressures of ethanol & cyclohexane at
# different temps using thermo parameters
# Psat = exp(A + B / T + C * log(T) + D * T ** E)
PsatEthanol = m.Intermediate(
m.exp(
PsatC2H5OH[0] +
PsatC2H5OH[1] / Temp +
PsatC2H5OH[2] * m.log(Temp) +
PsatC2H5OH[3] * Temp ** PsatC2H5OH[4]
)
)
PsatCyclohexane = m.Intermediate(
m.exp(
PsatC6H12[0] +
PsatC6H12[1] / Temp +
PsatC6H12[2] * m.log(Temp) +
PsatC6H12[3] * Temp ** PsatC6H12[4]
)
)
# Calculate Liquid Molar volumes of components 1 and 2:
# ethanol & cyclohexane, respectively
LMV_C2H5OH = 5.8492e-2
LMV_C6H12 = 0.10882
# Construct parameters for Wilson Model
L12 = m.Intermediate(LMV_C6H12 / LMV_C2H5OH \
* m.exp(-a12_ / (R * Temp)))
L21 = m.Intermediate(LMV_C2H5OH / LMV_C6H12 \
* m.exp(-a21_ / (R * Temp)))
L11 = 1.0
L22 = 1.0
x2 = m.Intermediate(1.0 - x1)
A1 = m.Intermediate(x1 + x2 * L12)
B1 = m.Intermediate(x1 * L11 / (x1 * L11 + x2 * L12) \
+ x2 * L21 / (x1 * L21 + x2 * L22))
A2 = m.Intermediate(x2 + x1 * L21)
B2 = m.Intermediate(x2 * L22 / (x2 * L22 + x1 * L21) \
+ x1 * L12 / (x2 * L12 + x1 * L11))
# Find Gammas for each component
G_C2H5OH = m.Intermediate(m.exp(1.0 - m.log(A1) - B1))
G_C6H12 = m.Intermediate(m.exp(1.0 - m.log(A2) - B2))
m.Equation(y2 == (G_C6H12 * PsatCyclohexane / P) * x2)
m.Equation(y1 == (G_C2H5OH * PsatEthanol / P) * x1)
# Options
m.options.IMODE = 2
m.options.EV_TYPE = 2 # Objective type, SSE
m.options.NODES = 3 # Collocation nodes
m.options.SOLVER = 3 # IPOPT
m.solve()
regressed_a12 = a12_.value[-1]
regressed_a21 = a21_.value[-1]
a12_Hysys = 2026.3
a21_Hysys = 390.4
x1_linspace = np.linspace(0, 1.0, 100)
y1_predicted = []
T_predicted = []
print("Activity Coefficient Parameters")
print('Model a12 a21')
print('Gekko: ' + str(round(regressed_a12,2)) \
+ str(round(regressed_a21,2)))
print('Hysys: ' + str(round(a12_Hysys,2)) \
+ str(round(a21_Hysys,2)))
def getValues(y1_and_Temp: list, a12: float, \
a21: float, P: float, x1: float):
y1 = y1_and_Temp[0]
T = y1_and_Temp[1]
y2 = 1 - y1
x2 = 1 - x1
A = y1*P - x1*Gamma(a12,a21,x1,T,1) \
* vaporP(T,PsatC2H5OH)
B = y2*P - x2*Gamma(a12,a21,x1,T,2) \
* vaporP(T,PsatC6H12)
return [A,B]
TempInterp = interp1d(x1_data, temp_data)
for i in range(len(x1_linspace)):
answers = fsolve(getValues, [x1_linspace[i],\
TempInterp(x1_linspace[i])],\
args = (regressed_a12, \
regressed_a21, \
P, x1_linspace[i]))
y1_predicted.append(answers[0])
T_predicted.append(answers[1])
GammaPred = Gamma(regressed_a12, regressed_a21, \
x1_data, temp_data, 1)
GammaPredHysys = Gamma(a12_Hysys, a21_Hysys, \
x1_data, temp_data, 2)
VP_Ethanol = vaporP(temp_data, PsatC2H5OH)
y1Predicted = (GammaPred * VP_Ethanol / P) * x1_data
y1PredictedHysys = (GammaPredHysys * VP_Ethanol / P) * x1_data
plt.figure(1)
plt.plot(x1_data,y1_data,"bx",label = "Measured Data")
plt.plot(x1_linspace,x1_linspace,"-",\
color="#8bd8bd", label="Reference",markersize=7)
plt.plot(x1_linspace,y1_predicted,"--",\
color="#243665",label = "Model",markersize=7)
plt.plot(x1_data, y1Predicted, "o",
color="#ec8b5e",label="Predicted",markersize=5)
plt.legend(loc = "best")
plt.xlabel("x")
plt.ylabel("y")
plt.grid()
plt.figure(2)
plt.plot(x1_data, temp_data, "go")
plt.plot(y1_data, temp_data, "bo")
plt.plot(x1_linspace, T_predicted, "r--",)
plt.plot(y1_predicted, T_predicted, "k--",)
plt.legend(["$x_1$ Measured", "$y_1$ Measured", \
"$x_1$ Predicted", "$y_1$ Predicted"])
plt.xlabel("$x, y$")
plt.ylabel("$Temp$")
plt.grid()
plt.show()