Importar las librerías

Cargar la tabla de datos


In [1]:
import scipy as sp
from scipy import optimize
from scipy.optimize import fsolve
import numpy as np
from matplotlib import pyplot
%matplotlib inline
import pandas as pd
from numpy import linalg as LA
from IPython.html import widgets
from IPython.display import display
from IPython.display import clear_output

# encoding: utf-8

from pandas import read_csv


/home/andres-python/anaconda3/lib/python3.5/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated since IPython 4.0. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
  "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)

Cargar la tabla de datos


In [2]:
f = pd.read_excel("PureFull.xls")
f.head()
data2 = pd.DataFrame(f)
data2 = data2.set_index('Name')
data2 = data2.ix[:, 1:12]
Etiquetas = data2.index.get_values()
Etiquetas


Out[2]:
array(['METHANE', 'ETHANE', 'PROPANE', ..., 'TITANIUM DIOXIDE',
       'PHOSPHORUS PENTOXIDE', 'PHOSPHORUS PENTOXIDE'], dtype=object)

In [3]:
Componentes_1 = widgets.SelectMultiple(
            description="Component 1",
            options=list(Etiquetas))
display(Componentes_1)

In [4]:
class Thermophysical_Properties():
    
    def __init__(self, nameData):
        self.nameData = nameData
        
    
    def cargar_Datos(self):
        f = pd.read_excel(self.nameData)
        f.head()
        data2 = pd.DataFrame(f)
        data2 = data2.set_index('Name')
        data2 = data2.ix[:, 1:12]
        self.Etiquetas = data2.index.get_values()
        
        print("Los datos del archivo: {0}, se han cargado correctamente !!!".format(self.nameData))
        
        return self.Etiquetas
    
    def seleccionar_Datos(self):
        Componentes_1 = widgets.SelectMultiple(
            description="Component 1",
            options=list(Etiquetas))
        display(Componentes_1)
        
        
    
    def mostrar_Datos(self):
        print ("Nombre componente: {0}".format(self.Etiquetas))
        


        
    
    def agregar_Datos(self):
        pass
    
    def borrar_Datos(self):
        pass
    
    def modificar_Datos(self):
        pass
    
    def crear_Datos(self):
        pass

In [5]:
nameData = "PureFull.xls"

propiedades = Thermophysical_Properties(nameData)
propiedades.cargar_Datos()
propiedades.mostrar_Datos()
propiedades.seleccionar_Datos()


Los datos del archivo: PureFull.xls, se han cargado correctamente !!!
Nombre componente: ['METHANE' 'ETHANE' 'PROPANE' ..., 'TITANIUM DIOXIDE'
 'PHOSPHORUS PENTOXIDE' 'PHOSPHORUS PENTOXIDE']

In [6]:
class System_Conditions():
    
    def __init__(self, Temperature, Pressure, Volume, Mole_fraction, Model_fluid, Model_solid):
        self.Temperature = Temperature
        self.Mole_fraction = Mole_fraction
        pass
    
    def normalizar(self):
        self.Mole_fraction_normal = Mole_fraction / sum(self.Mole_fraction)
        return self.Mole_fraction_normal
    
    def convertir(self):
        pass
        
        
    
class Componentes(Thermophysical_Properties, System_Conditions):
    """
    Las variables aux_ se utilizan para presentar de forma más clara y acotada
    las expresiones necesarias en los calculos. Estas, se numeran de acuerdo al orden de
    aparición dentro de una clase.
    
    """
    
    def __init__(self):
        pass
    
    def cal_SRK_model(self):
        # Soave-Redlich-Kwong (SRK)
        self.s1, self.s2 = 1, 2
        self.m = 0.480 + 1.574 * self.w - 0.175 * self.w ** 2
        self.ac = 0.077796070 * self.R ** 2, self.Tc ** 2 / self.Pc
        self.bc = 0.086640 * self.R * self.Tc / self.Pc        
        
        return self.m, self.ac, self.bc
    
    def cal_PR_model(self):
        # Peng-Robinson (PR)
        self.s1, self.s2 = 1 + 2 ** 0.5, 1 - (2 ** 0.5)
        self.m = 0.37464 + 1.54226 * self.w - 0.26992 * self.w ** 2
        self.ac = 0.45723553 * self.R ** 2 * self.Tc ** 2 / self.Pc
        self.bc = 0.077796070 * self.R * self.Tc / self.Pc            
        
        self.alfa = (1 + self.m * (1 - (self.T / self.Tc) ** 0.5)) ** 2
        aux_1 = - (self.m / self.T) * (self.T / self.Tc) ** 0.5
        aux_2 = (self.m * (- (self.T / self.Tc) ** 0.5 + 1) + 1)
        self.dalfadT = aux_1 * aux_2
        
        aux_3 = 0.5 * self.m ** 2 * (self.T / self.Tc) ** 1.0 / self.T ** 2
        aux_4 = (self.m * (- (self.T / self.Tc) ** 0.5 + 1) + 1) / self.T ** 2
        aux_5 = 0.5 * self.m * (self.T / self.Tc) ** 0.5 * aux_4
        
        self.d2alfaT2 = aux_3 + aux_5
        
        self.a_ii = self.ac * self.alfa
        self.b_ii = self.bc
        self.da_iidT = self.ac * self.dalfadT
        d2adT2_puros = self.ac * self.d2alfaT2
        
        return self.m, self.a_ii, self.b_ii
    
    def cal_RKPR_model(self):
        pass
    
    def build_component(self):
        
        if self.eq == "SRK":
            # Soave-Redlich-Kwong (SRK)
            self.component = self.cal_SRK_model()
        elif self.eq == "PR":
            # Peng-Robinson (PR)
            self.component = self.cal_PR_model()
        elif self.eq == "RKPR":
            # (RKPR)
            #self.component = self.cal_RKPR_model()
            print ("No actualizada, intentalo de nuevo !!! ")
        else:
            print ("Che boludo... Modelo no valido, intentalo de nuevo !!! ")

In [7]:
Componentes_1 = widgets.SelectMultiple(
    description="Component 1",
    options=list(Etiquetas))

Componentes_2 = widgets.SelectMultiple(
    description="Component 2",
    options=list(Etiquetas))

button = widgets.Button(description="Upload Data")

def cargarDatos(b):
    clear_output()
    print("Component 1: ", Componentes_1.value)
    Nombre = Componentes_1.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_1 = Propiedades[0]
    Temperatura_Critica_1 = Propiedades[1]
    Presion_Critica_1 = Propiedades[2]
    Z_Critico_1 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor = ", Factor_Acentrico_1)
    print ("Critical Temperature = ", Temperatura_Critica_1, "K")
    print ("Critical Pressure = ", Presion_Critica_1, "bar")
    print ("Z_Critical = ", Z_Critico_1, "\n")
    
    
    print("Component 2: ", Componentes_2.value)
    Nombre = Componentes_2.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_2 = Propiedades[0]
    Temperatura_Critica_2 = Propiedades[1]
    Presion_Critica_2 = Propiedades[2]
    Z_Critico_2 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor  = ", Factor_Acentrico_2)
    print ("Critical Temperature = ", Temperatura_Critica_2, "K")
    print ("Critical Pressure = ", Presion_Critica_2, "bar")
    print ("Z_Critical = ", Z_Critico_2)
    
    global TcDato, PcDato, wDato
    
    TcDato = np.array([Temperatura_Critica_1, Temperatura_Critica_2])
    PcDato = np.array([Presion_Critica_1, Presion_Critica_2])
    wDato = np.array([Factor_Acentrico_1, Factor_Acentrico_2])


button.on_click(cargarDatos)

In [8]:
class Parameters_BD():
    
    def __init__(self):
        pass
    
    def cal_parameters_ij(self):           

        if self.nC > 1:
            self.aij = np.ones((len(self.ni), len(self.ni)))
            self.bij = np.ones((len(self.ni), len(self.ni)))
            self.daijdT = np.ones((len(self.ni), len(self.ni)))

            for j in range(self.nC):
                for i in range(self.nC):
                    self.aij[i, j] = (self.a_ii[i] * self.a_ii[j]) ** 0.5
                    self.bij[i, j] = (self.b_ii[i] + self.b_ii[j]) / 2
                    self.bij[i, j] = self.bij[i, j]
                    self.daijdT[i, j] = (self.da_iidT[i] * self.da_iidT[j]) ** 0.5

            for i in range(self.nC):
                for  j in range(self.nC):
                    if i == j:
                        self.aij[i, j] = self.a_ii[i] * (1 - self.kij[i, j])
                        self.daijdT[i, j] = self.da_iidT[i] * (1 - self.kij[i, j])
                    elif i != j:
                        self.aij[i, j] = self.aij[i, j] * (1 - self.kij[i, j])
                        self.daijdT[i, j] = self.daijdT[i, j] * (1 - self.kij[i, j])
                   
        if self.nC == 1:
            return self.a_ii, self.b_ii, self.da_iidT
        else:
            return self.aij, self.bij, self.daijdT

    def cal_parameter_D(self):
        if self.nC == 1:
            self.D = self.ni ** 2 * self.a_ii
            self.Di = 2 * self.ni * self.a_ii
        else:
            di = np.ones((len(self.ni), len(self.ni)))
            self.Di = np.ones((len(self.ni)))
            self.D = np.ones((len(self.ni)))
            for i in range(self.nC):
                for j in range(self.nC):
                    di[i, j] = self.ni[j] * self.aij[i, j]
                    self.Di[i] = 2 * np.sum(di[i, :])
            self.D = 0.5 * np.sum(self.ni * self.Di)

        return self.D
    
    def cal_parameter_delta_1(self):
        
        if self.nC == 1:
            self.D1m = np.zeros((len(self.ni)-1))
            self.dD1i = np.ones((len(self.ni)))
            self.dD1ij = np.ones((len(self.ni), len(self.ni)))
            
            for i in range(self.nC):
                self.D1m = self.D1m + self.ni[i] * self.delta_1[i]
            
            self.D1m = self.D1m / self.nT
            
        else:
            self.D1m = np.zeros((len(self.ni)-1))
            self.dD1i = np.ones((len(self.ni)))
            self.dD1ij = np.ones((len(self.ni), len(self.ni)))
            
            for i in range(self.nC):
                self.D1m = self.D1m + self.ni[i] * self.delta_1[i]
            
            self.D1m = self.D1m / self.nT
            
            for i in range(self.nC):
                self.dD1i[i] = (self.delta_1[i] - self.D1m) / self.nT
                for j in range(self.nC):
                    self.dD1ij[i,j] = (2.0 * self.D1m - self.delta_1[i] - self.delta_1[j]) / self.nT ** 2
                    
        return self.D1m, self.dD1i, self.dD1ij
                    
    def cal_parameter_B(self):
        if self.nC == 1:
            self.B = self.ni * self.b_ii
        else:
            self.aux = np.zeros((len(self.ni)))
            for i in range(self.nC):
                for j in range(self.nC):
                    self.aux[i] = self.aux[i] + self.ni[j] * self.bij[i, j]

            self.B = np.sum(self.ni * self.b_ii)
            #print("B = ", self.B)

        return self.B

In [9]:
class Fugacidad():

    def __init__(self, eq, w, Tc, Pc, Tr, R, ep, ni, nT, nC, V, T, P, kij, lij, delta_1, k, Avsl):
        self.eq = eq
        self.w = w
        self.Tc = Tc
        self.Pc = Pc
        self.Tr = Tr
        self.R = R
        self.ep = ep
        self.ni = ni
        self.nT = nT
        self.nC = nC
        self.V = V
        self.T = T
        self.P = P
        self.kij = kij
        self.lij = lij
        self.delta_1 = delta_1
        self.k = k
        self.Avsl = Avsl
        
        if self.eq == "SRK":
            # Soave-Redlich-Kwong (SRK)
            self.s1, self.s2 = 1, 2
            self.m = 0.480 + 1.574 * self.w - 0.175 * self.w ** 2
            self.ac = 0.077796070 * self.R ** 2, self.Tc ** 2 / self.Pc
            self.bc = 0.086640 * self.R * self.Tc / self.Pc
        elif self.eq == "PR":
            # Peng-Robinson (PR)
            self.s1, self.s2 = 1 + 2 ** 0.5, 1 - (2 ** 0.5)
            self.m = 0.37464 + 1.54226 * self.w - 0.26992 * self.w ** 2
            self.ac = 0.45723553 * self.R ** 2 * self.Tc ** 2 / self.Pc
            self.bc = 0.077796070 * self.R * self.Tc / self.Pc            
           
            self.alfa = (1 + self.m * (1 - (self.T / self.Tc) ** 0.5)) ** 2
            self.dalfadT = - (self.m / self.T) * (self.T / self.Tc) ** 0.5 * (self.m * (- (self.T / self.Tc) ** 0.5 + 1) + 1)
            ter_1 = 0.5 * self.m ** 2 * (self.T / self.Tc) ** 1.0 / self.T ** 2
            ter_2 = 0.5 * self.m * (self.T / self.Tc) ** 0.5 * (self.m * (- (self.T / self.Tc) ** 0.5 + 1) + 1) / self.T ** 2
            
            self.d2alfaT2 = ter_1 + ter_2
            self.a_ii = self.ac * self.alfa
            self.b_ii = self.bc
            
            self.da_iidT = self.ac * self.dalfadT
            d2adT2_puros = self.ac * self.d2alfaT2

        elif self.eq == "RKPR":
            # (RKPR)
            print ("No actualizada, intentalo de nuevo !!! ")            

        else:
            print ("Che boludo... Modelo no valido, intentalo de nuevo !!! ")


    def parametros(self):           

        if self.nC > 1:
            self.aij = np.ones((len(self.ni), len(self.ni)))
            self.bij = np.ones((len(self.ni), len(self.ni)))
            self.daijdT = np.ones((len(self.ni), len(self.ni)))

            for j in range(self.nC):
                for i in range(self.nC):
                    self.aij[i, j] = (self.a_ii[i] * self.a_ii[j]) ** 0.5
                    self.bij[i, j] = (self.b_ii[i] + self.b_ii[j]) / 2
                    self.bij[i, j] = self.bij[i, j]
                    self.daijdT[i, j] = (self.da_iidT[i] * self.da_iidT[j]) ** 0.5

            for i in range(self.nC):
                for  j in range(self.nC):
                    if i == j:
                        self.aij[i, j] = self.a_ii[i] * (1 - self.kij[i, j])
                        self.daijdT[i, j] = self.da_iidT[i] * (1 - self.kij[i, j])
                    elif i != j:
                        self.aij[i, j] = self.aij[i, j] * (1 - self.kij[i, j])
                        self.daijdT[i, j] = self.daijdT[i, j] * (1 - self.kij[i, j])
                   
        if self.nC == 1:
            return self.a_ii, self.b_ii, self.da_iidT
        else:
            return self.aij, self.bij, self.daijdT

    def parametro_D(self):
        if self.nC == 1:
            self.D = self.ni ** 2 * self.a_ii
            self.Di = 2 * self.ni * self.a_ii
        else:
            di = np.ones((len(self.ni), len(self.ni)))
            self.Di = np.ones((len(self.ni)))
            self.D = np.ones((len(self.ni)))
            for i in range(self.nC):
                for j in range(self.nC):
                    di[i, j] = self.ni[j] * self.aij[i, j]
                    self.Di[i] = 2 * np.sum(di[i, :])
            self.D = 0.5 * np.sum(self.ni * self.Di)

        return self.D
    
    def parametro_delta_1(self):
        
        if self.nC == 1:
            self.D1m = np.zeros((len(self.ni)-1))
            self.dD1i = np.ones((len(self.ni)))
            self.dD1ij = np.ones((len(self.ni), len(self.ni)))
            
            for i in range(self.nC):
                self.D1m = self.D1m + self.ni[i] * self.delta_1[i]
            
            self.D1m = self.D1m / self.nT
            
        else:
            self.D1m = np.zeros((len(self.ni)-1))
            self.dD1i = np.ones((len(self.ni)))
            self.dD1ij = np.ones((len(self.ni), len(self.ni)))
            
            for i in range(self.nC):
                self.D1m = self.D1m + self.ni[i] * self.delta_1[i]
            
            self.D1m = self.D1m / self.nT
            
            for i in range(self.nC):
                self.dD1i[i] = (self.delta_1[i] - self.D1m) / self.nT
                for j in range(self.nC):
                    self.dD1ij[i,j] = (2.0 * self.D1m - self.delta_1[i] - self.delta_1[j]) / self.nT ** 2
                    
        return self.D1m, self.dD1i, self.dD1ij
                    
    def parametro_B(self):
        if self.nC == 1:
            self.B = self.ni * self.b_ii
        else:
            self.aux = np.zeros((len(self.ni)))
            for i in range(self.nC):
                for j in range(self.nC):
                    self.aux[i] = self.aux[i] + self.ni[j] * self.bij[i, j]

            self.B = np.sum(self.ni * self.b_ii)
            #print("B = ", self.B)

        return self.B

    def presion(self):
        '''
        Con el metodo presion(), se calcula la Presión P(T, V, N) del sistema
        para una temperatura T, cantidad de moles N y un volumen V
        R = Constante universal de los gases
        nT = Número total de moles en el sistema
        Pcal = Peos = Presión calculada con la ecuación de estado
        Arv = Primera derivada parcial de la energía de Helmholz con respecto al
        volumen V, a T y N constantes
        '''
        self.gv = self.R * self.B / (self.V * (self.V - self.B))
        self.fv = - 1 / ((self.V + self.s1 * self.B) * (self.V + self.s2 * self.B))
        self.ArV = -self.nT * self.gv * self.T - self.D * self.fv
        self.Pcal = self.nT * self.R * self.T / self.V - self.ArV   

        return self.Pcal

    def dP_dV(self):
        self.dPdV = -self.ArV2 - self.R * self.T * self.nT / self.V ** 2
        return self.dPdV

    def Z_factor(self):
        self.Z = (self.P * self.V) / (self.nT * self.R * self.T)
        return self.Z

    def P_ideal(self):
        self.Pxi = (self.ni * self.P) / self.nT
        return self.Pxi

    def dF_dV(self):
        '''
        Primera derivada de F con respecto al volumen Ecu. (68)
        '''
        self.gv = self.R * self.B / (self.V * (self.V - self.B))
        self.fv = - 1 / ((self.V + self.s1 * self.B) * (self.V + self.s2 * self.B))
        self.ArV = -self.nT * self.gv * self.T - self.D * self.fv
        return self.ArV

    def dF_dVV(self):
        '''
        Segunda derivada de F con respecto al volumen Ecu. (74)
        '''
        self.gv2 = self.R * (1 / self.V ** 2 - 1 / (self.V - self.B) ** 2)
        self.fv2 = (- 1 / (self.V + self.s1 * self.B) ** 2 + 1 / (self.V + self.s2 * self.B) ** 2) / self.B / (self.s1 - self.s2)
        self.ArV2 = - self.nT * self.gv2 * self.T - self.D * self.fv2
        return self.ArV2

    def volumen_1(self):
        '''
        Calculo del volumen V(T,P,n) del fluido a una temperatura T, presión P
        y número de moles totales nT especificados.
        Se utiliza el método de Newton con derivada de la función analitica.
        Pendiente cambiar por una función de Scipy.
        '''
        self.V = 1.05 * self.B          
        lnP = np.log(self.P)
        Pite = self.presion()
        lnPcal = np.log(Pite)
        h = lnP - lnPcal
        errorEq = abs(h)
        i = 0
        s = 1.0

        while errorEq > self.ep:
            self.parametro_D()
            self.parametro_B()
            self.dF_dV()
            self.dF_dVV()
            dPite = self.dP_dV()
            Pite = self.presion()
            lnPcal = np.log(Pite)
            h = lnP - lnPcal
            dh = -dPite
            self.V = self.V - s * h / dh
            errorEq = abs(h)
            i += 1
            if i >= 900:
                pass
                #break

        return self.V

    def funcion_energia_F(self):
        self.g = self.R * np.log(1 - self.B / self.V)
        self.bv = self.B / self.V
        self.f = np.log((self.V + self.s1 * self.B) / (self.V + self.s2 * self.B)) / self.B / (self.s1 - self.s2)
        self.Ar = -self.nT * self.g * self.T - self.D * self.f
        return self.g, self.f, self.Ar, self.bv
    
    def tomar_B(self):
        print ("tomando B =", self.B)
        return self.B + 10
    
    def derivadas_delta_1(self):
        auxD2 = (1 + 2 / (1 + self.s1) ** 2)
        
        como_1 = (1 / (self.V + self.s1 * self.B) + 2 / (self.V + self.s2 * self.B) / (1 + self.s1) ** 2)
        como_2 = self.f * auxD2
        self.fD1 = como_1 - como_2
        self.fD1 = self.fD1/(self.s1 - self.s2)
        
        return self.fD1

    def primeras_derivadas1(self):

        if self.nC == 1:
            AUX = self.R * self.T / (self.V - self.B)
            self.fB = -(self.f + self.V * self.fv) / self.B
            self.FFB = self.nT * AUX - self.D * self.fB
            self.Di = 2 * self.nT * self.ac * self.alfa
            self.Bi = self.bc
            
            if self.eq != "RKPR":
                self.Arn = -self.g * self.T + self.FFB * self.Bi - self.f * self.Di
            else:
                self.Arn = -self.g * self.T + self.FFB * self.Bi - self.f * self.Di - self.D * self.fD1 * self.dD1i
        else:
            # Derivando la ecuación (64) se obtiene la ecuación eq (106)
            self.Bi = np.ones((len(self.ni)))
            for i in range(self.nC):
                self.Bi[i] = (2 * self.aux[i] - self.B) / self.nT

            AUX = self.R * self.T / (self.V - self.B)
            self.fB = -(self.f + self.V * self.fv) / self.B
            self.FFB = self.nT * AUX - self.D * self.fB
            
            if self.eq != "RKPR":
                self.Arn = -self.g * self.T + self.FFB * self.Bi - self.f * self.Di
            else:
                auxD2 = (1 + 2 / (1 + self.s1) ** 2)
                print("B delta1 = ", self.B)
                co_1 = (1 / (self.V + self.s1 * self.B) + 2 / (self.V + self.s2 * self.B) / (1 + self.s1) ** 2)
                co_2 = self.f * auxD2
                self.fD1 = co_1 - co_2
                self.fD1 = self.fD1/(self.s1 - self.s2)
                self.Arn = -self.g * self.T + self.FFB * self.Bi - self.f * self.Di - self.D * self.fD1 * self.dD1i

        return self.Arn, self.Arn, self.Arn
    

    def coeficientes_fugacidad(self):
        self.Z = self.Z_factor()
        self.lnOi = self.Arn / (self.R * self.T) - np.log(self.Z)
        self.Oi = np.exp(self.lnOi)
        return self.Oi

    def fugacidad(self):
        self.Z = self.Z_factor()
        self.Pxi = self.P_ideal()
        self.lnFi = self.Arn / (self.R * self.T) - np.log(self.Z) + np.log(self.Pxi)
        self.Fi = np.exp(self.lnFi)
        self.PHILOG = self.Arn / (self.R * self.T) - np.log(self.Z)
        self.PHILOG_i = self.Arn - np.log(self.Z)
        self.FUGLOG = self.Arn / (self.R * self.T) + np.log(self.ni) + np.log((self.nT * self.R * self.T) / self.V)
        return self.Fi

    def exp_sol(self):
        '''
        Este método calcula el factor de corrección de la fugacidad del
        componente fluido para determinar la fugacidad del mismo componente
        en estado sólido.
        Fugacidad del sólido puro
        fi_s(T, P) = fi_l(T, P) * EXP(T, P)
        '''
        Tfus = 323.75
        # Temperatura de fusion de n-tetracosane
        # Unidad de Ti_f en Kelvin
        par_sol = np.array([[-176120.0, 8196.20, -55.911, 0.19357, -0.0002235],
                            [-1.66e6, 8.31e3, 0.0, 0.0, 0.0]])
        par_liq = np.array([[423160.0, 1091.9, 0.0, 0.0, 0.0],
                            [7.01e5, 1.47e3, 0.0, 0.0, 0.0]])
        #print ("par_sol", par_sol)
        #print ("par_liq", par_liq)
        # Las unidades de Cp están en J/Kmol.K
        Cp_solido = par_sol[:, 0] + par_sol[:, 1] * T + par_sol[:, 2] * T ** 2 + par_sol[:, 3] * T ** 3 + par_sol[:, 4] * T ** 4
        #print ("Cp_solido", Cp_solido)
        Cp_liquido= par_liq[:, 0] + par_liq[:, 1] * T + par_liq[:, 2] * T ** 2 + par_liq[:, 3] * T ** 3 + par_liq[:, 4] * T ** 4
        #print ("Cp_liquido", Cp_liquido)
        DeltaCp = (Cp_solido - Cp_liquido) * (1.0 / 1000)
        print ("Delta Cp", DeltaCp)

        #Unidades de Delta H de fusión en Kcal/mol
        DeltaH_f = np.array([13.12, 21.23]) * (1000 / 1.0) * (4.18 / 1.0)
        #print ("Delta H de fusion", DeltaH_f)
        T_f = np.array([323.75, 349.05])
        #print ("Temperaturas de fusion = ", T_f)

        Rp = 8.314
        A = (DeltaH_f / (Rp * Tfus)) * (1 - (Tfus / T))
        B = (DeltaCp / Rp) * (1 - (Tfus / T))
        C = (DeltaCp / Rp) * np.log(Tfus / T)
        self.EXP = np.exp(A - B - C)

        print ("A = ", A)
        print ("B = ", B)
        print ("C = ", C)
        print ("EXP = ", self.EXP)

        return self.EXP
    
    def exp_sol_1(self):
        '''
        Este método calcula el factor de corrección de la fugacidad del
        componente fluido para determinar la fugacidad del mismo componente
        en estado sólido.
        Fugacidad del sólido puro
        fi_s(T, P) = fi_l(T, P) * EXP(T, P)
        '''
        Tpt = 323.75
        Ppt = 1.38507E-8
        R = 8.314472
        AH = 54894000
        Av = -0.0376300841 #m3/kmol
        
        a = ((AH / (R * Tpt)) * (1 - (Tpt / self.T))) / 1000
        b = ((Av / (R * self.T)) * (self.P - Ppt)) * 100
        self.EXP_1 = a + b
        
        return self.EXP_1
    
    def exp_sol_3(self):
        '''
        Este método calcula el factor de corrección de la fugacidad del
        componente fluido para determinar la fugacidad del mismo componente
        en estado sólido.
        Fugacidad del sólido puro
        fi_s(T, P) = fi_l(T, P) * EXP(T, P)
        '''
        # [=] K
        # [=] bar
        # [m3 / Kmol]
        # Constante R [=] 0.08314472 bar.l/(mol.K)
        
        Tpt = 323.75
        Ppt = 3.2015002E-8
        #self.Avsl = -0.0565500835
        
        c1 = -14213.5004
        c2 = 605153.4382
        c3 = -591592.556
        
        R = 0.08314472
        
        A1 = c1 * (1 - Tpt / self.T)
        A2 = c2 * (-1 + Tpt / self.T + np.log(self.T / Tpt))
        A3 = c3 * (-1 + self.T / (2 * Tpt) + Tpt / (2 * self.T)) + (Tpt / self.T) * (self.P - Ppt)
        
        FE = (self.Avsl / (self.R * self.T)) * (A1 + A2 + A3)
        self.EXP_3 = np.exp(FE)    
        return self.EXP_3
        

    def fluido(self):
        ab = self.parametros()
        D = self.parametro_D()
        B = self.parametro_B()
        Vol_1 = self.volumen_1()
        F = self.funcion_energia_F()
        dF = self.primeras_derivadas1()
        Z = self.Z_factor()
        Zcal = (self.P * Vol_1) / (self.nT * self.R * self.T)
        Pq = self.presion()
        self.Fug = self.fugacidad()
        self.CoeFug = self.coeficientes_fugacidad()
        return self.Fug

    def solido(self):
        if self.nC == 1:
            Fug = self.fluido()
            #EXP = self.exp_sol()
            #EXP = self.exp_sol_1()
            EXP = self.exp_sol_3()
            
            FugS = Fug[0] * EXP
        else:
            print ("Aún no se qué hacer para una mezcla de sólidos !!!")
            FugS = 1

        return FugS
    
#----------------
def calculaFugacidad(x, Pe, nif, nCf, eq, TcDato, PcDato, wDAto, Avsl):
    #---------------------------------------------------------------------------
    # Temperatura en [=] K
    # Presión en [=] bar
    # Constante R [=] 0.08314472 bar.l/(mol.K)
    # x = variable que se cálcula, puede ser T ó P para el equilibrio sólido-fluido
    # Pe = Presión del sistema especificada
    # nif = número de moles del componente (i) en cada fase (f)
    # nCf = número de componentes en una fase (f)
    # eq = modelo de ecuación de estado, SRK, PR, RKPR
    # TcDato = Temperatura critica de la "base de datos"
    # PcDato = Presión critica de la "base de datos"
    # wDato = Factor acentrico de la "base de datos"
    # Avsl = Delta de volumen sólido-fluido
   
    # ep = Criterio de convergencia del método def volumen_1(self, P)
    
    T = x # 335.42 # x # 366.78 # 356.429 # 335.42 # 348.89 #327.0
    #print("Temperatura = ", T)
    P = Pe # 2575.0 # 2064.7 # 1524.4 #1164.2 # 865.0 
    # 560.3 # x #1054.6 #1560.3 # 2064.7 # 1524.4 # 560.3 # 1164.2 #865.0
    R = 0.08314472
    ep = 1e-5#1e-6
    #---------------------------------------------------------------------------    
    Tcm = TcDato
    Pcm = PcDato
    wm = wDato
        
    nC = nCf
    
    if nC == 1:
        #print ("...............................................................")
        
        #ni = nif
        ni = np.array([1.0])
        
        #print ("Número de moles = ", ni)
        # C24
        kij = 0.0
        lij = 0.0
        
        # Metano - Etano
        delta_1 = np.array([0.85])        
        k = np.array([1.50758])
        #C24
        Tc = Tcm[1]
        Pc = Pcm[1]
        w = wm[1]     
        print ("...............................................................")
    elif nC == 2:
        # metano - C24
        #ni = np.array([1-nif, nif])
        ni = nif #np.array([1-nif, nif])
        
        #ni = np.array([1 - 0.901, 0.901])
        #---------------------------------
        
        
        #ni = np.array([1 - 0.26, 0.26])
        
        #ni = np.array([1 - 0.104, 0.104])
        #print ("Número de moles = ", ni)

        kij = np.array([[0.000000, 0.083860],
                        [0.083860, 0.000000]])
        
        kij = np.array([[0.000000, 0.059600],
                        [0.059600, 0.000000]])
        

        lij = 0.0132
        
        
        #kij = np.array([[0.000000, 0.00],
        #                [0.00, 0.000000]])
        
        #lij = 0.0
        
        
        # Metano - C24
        delta_1 = np.array([0.85, 2.40])        
        k = np.array([1.50758, 4.90224])
        
        # metano sigma1 = 0.9253, sigma = 0.85, k = 1.49345, k = 1.50758
        # C24 sigma = 2.40 k = 4.90224

        Tc = Tcm
        Pc = Pcm
        w = wm
        print ("Temperatura Critica = ", Tc, "K")
        print ("Presión Critica = ", Pc, "bar")
        print ("Factor Acentrico = ", w)
        #print ("...............................................................")

    # Tempertura reducidad
    Tr = T / Tc
    # C24 puro
    V = 0.141604834257319
    nT = np.sum(ni)

    fugacidad = Fugacidad(eq, w, Tc, Pc, Tr, R, ep, ni, nT, nC, V, T, P, kij, lij, delta_1, k, Avsl)
    
    print(fugacidad.exp_sol_3())
            
    if nC == 1:
        SOL = fugacidad.solido()
        
        return SOL
    else:
        flu_1 = fugacidad.fluido()
        return flu_1


#----------------

In [10]:
def equilibrioSF(x, Pe, nif, n1, n2, Avsl):
    
    # fugacidad del sólido puro
    FugS = calculaFugacidad(x, Pe, nif, n1, eq, TcDato, PcDato, wDato, Avsl)
    print(eq, TcDato, PcDato, wDato, Avsl)
    # fugacidad del fluido pesado en la mezcla fluida
    FugF = calculaFugacidad(x, Pe, nif, n2, eq, TcDato, PcDato, wDato, Avsl)
    
    # Función de igualdad de fugacidades del sólido y el fluido
    eqSF = np.abs(np.abs(np.log(FugS)) - np.abs(np.log(FugF[1])))
    print ("-"*80)
    print ("ln(Fugacidad Sólido) = ", np.log(FugS))
    print ("ln(Fugacidad Fluido) = ", np.log(FugF[1]))
    print ("ln(Fugacidad Sólido) - ln(Fugacidad Fluido) = ", eqSF)   
    
    return eqSF


eq = 'PR'
Avsl = -0.0565500835
#Avsl = -0.09605965500835

initial_temperature = [346.5] # T [=] K
initial_pressure = 136.9 # [=] bar

Tcal = fsolve(equilibrioSF,initial_temperature,args=(initial_pressure, 1, 2, Avsl), xtol=1e-4)
print(Tcal, "K")


---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-10-6cc4466ce984> in <module>()
     24 initial_pressure = 136.9 # [=] bar
     25 
---> 26 Tcal = fsolve(equilibrioSF,initial_temperature,args=(initial_pressure, 1, 2, Avsl), xtol=1e-4)
     27 print(Tcal, "K")

/home/andres-python/anaconda3/lib/python3.5/site-packages/scipy/optimize/minpack.py in fsolve(func, x0, args, fprime, full_output, col_deriv, xtol, maxfev, band, epsfcn, factor, diag)
    144                'diag': diag}
    145 
--> 146     res = _root_hybr(func, x0, args, jac=fprime, **options)
    147     if full_output:
    148         x = res['x']

/home/andres-python/anaconda3/lib/python3.5/site-packages/scipy/optimize/minpack.py in _root_hybr(func, x0, args, jac, col_deriv, xtol, maxfev, band, eps, factor, diag, **unknown_options)
    210     if not isinstance(args, tuple):
    211         args = (args,)
--> 212     shape, dtype = _check_func('fsolve', 'func', func, x0, args, n, (n,))
    213     if epsfcn is None:
    214         epsfcn = finfo(dtype).eps

/home/andres-python/anaconda3/lib/python3.5/site-packages/scipy/optimize/minpack.py in _check_func(checker, argname, thefunc, x0, args, numinputs, output_shape)
     24 def _check_func(checker, argname, thefunc, x0, args, numinputs,
     25                 output_shape=None):
---> 26     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
     27     if (output_shape is not None) and (shape(res) != output_shape):
     28         if (output_shape[0] != 1):

TypeError: equilibrioSF() missing 1 required positional argument: 'Avsl'

In [11]:
def equilibrioSF(x, Pe, nif, n1, n2, Avsl):
    
    # fugacidad del sólido puro
    FugS = calculaFugacidad(x, Pe, nif, n1, eq, TcDato, PcDato, wDato, Avsl)
    print(eq, TcDato, PcDato, wDato, Avsl)
    # fugacidad del fluido pesado en la mezcla fluida
    FugF = calculaFugacidad(x, Pe, nif, n2, eq, TcDato, PcDato, wDato, Avsl)
    
    # Función de igualdad de fugacidades del sólido y el fluido
    eqSF = np.abs(np.abs(np.log(FugS)) - np.abs(np.log(FugF[1])))
    print ("-"*80)
    print ("ln(Fugacidad Sólido) = ", np.log(FugS))
    print ("ln(Fugacidad Fluido) = ", np.log(FugF[1]))
    print ("ln(Fugacidad Sólido) - ln(Fugacidad Fluido) = ", eqSF)   
    
    return eqSF


eq = 'PR'
#Avsl = -0.0565500835
#Avsl = -0.09605965500835

#initial_temperature = [346.5] # T [=] K
#initial_pressure = 136.9 # [=] bar

#Tcal = fsolve(equilibrioSF,initial_temperature,args=(initial_pressure, 1, 2, Avsl), xtol=1e-4)
#print(Tcal, "K")

t_exp = [323.65, 326.04, 326.43, 328.12, 329.45, 329.89, 333.43, 335.12, 340.19, 344.58, 346.65, 352.53, 362.45, 362.76, 371.82, 379.74]
temp = np.array(t_exp)

p_exp = [1, 101.0, 136.9, 183.8, 266.2, 266.8, 426.9, 480.3, 718.9, 912.5, 1010.6, 1277.8, 1778.0, 1825.1, 2323.4, 2736.1]
pres= np.array(p_exp)

pos = np.arange(len(pres))
Tcal = np.ones((len(pres)))
Tcal

Tres = np.array([ 322.65861561,  324.91946742,  325.73456905,  326.80151121,
        328.68045402,  328.69415114,  332.3526483 ,  333.57248076,
        338.99640222,  343.33723415,  345.50684642,  351.28742799,
        361.49784425,  362.4145721 ,  371.63445321,  378.63493779])

Tcal - temp 

Avsl = -0.32595074
Avsl


class Flash():

    def __init__(self, zi_F, temperature_f, pressure_f, TcDato_f, PcDato_f, wDato_f):
        self.zi = zi_F
        self.T = temperature_f
        self.P = pressure_f
        self.Tc = TcDato_f
        self.Pc = PcDato_f
        self.w = wDato_f        
        
    def wilson(self):
        # Ecuación wilson
        lnKi = np.log(self.Pc / self.P) + 5.373 * (1 + self.w) * (1 - self.Tc / self.T)
        self.Ki = np.exp(lnKi)
        return self.Ki

    def beta(self):
        # Estimación de la fracción de fase de vapor en el sistema
        self.Ki = self.wilson()
        #Bmin = np.divide((self.Ki * self.zi - 1), (self.Ki - 1))
        Bmin = (self.Ki * self.zi - 1) / (self.Ki - 1)
        
        #print (("Bmin_inter = ", Bmin))
        
        Bmax = (1 - self.zi) / (1 - self.Ki)
        #print (("Bmax_inter = ", Bmax))
        self.Bini = (np.max(Bmin) + np.min(Bmax)) / 2
        print("inib =", self.Bini)
        return self.Bini

    def rice(self):
        # Ecuación de Rachford-Rice para el equilibrio líqudo-vapor
        self.fg = np.sum(self.zi * (self.Ki - 1) / (1 - self.Bini + self.Bini * self.Ki))
        self.dfg = - np.sum(self.zi * (self.Ki - 1) ** 2 / (1 - self.Bini + self.Bini * self.Ki) ** 2)
        #print g, dg
        return self.fg, self.dfg
    
    def composicion_xy(self):
        # Ecuación de Rachford-Rice para calcular la composición del equilibrio líqudo-vapor
        self.xi = self.zi / (1 - self.Bini + self.Bini * self.Ki)
        self.yi = (self.zi * self.Ki) / (1 - self.Bini + self.Bini * self.Ki)
        self.li = (self.zi * (1 - self.Bini)) / (1 - self.Bini + self.Bini * self.Ki)
        self.vi = (self.zi * self.Bini * self.Ki) / (1 - self.Bini + self.Bini * self.Ki)

        return self.xi, self.yi, self.li, self.vi

    def flash_ideal(self):
        # Solución del flash (T,P,ni) isotermico para Ki_(T,P)
        self.Bini = self.beta()
        self.Ki = self.wilson()
        # print ("Ki_(P, T) = ", self.Ki)
        Eg = self.rice()
        errorEq = abs(Eg[0])
        # Especificaciones del método Newton precario, mientras se cambia por una librería Scipy
        i, s, ep = 0, 1, 1e-5

        while errorEq > ep:
            Eg = self.rice()
            self.Bini = self.Bini - s * Eg[0] / Eg[1]
            errorEq = abs(Eg[0])
            i += 1
            if i >= 50:
                break

        xy = self.composicion_xy()        
        print ("-"*53, "\n", "-"*18, "Mole fraction", "-"*18, "\n","-"*53)
        print ("\n", "-"*13, "Zi phase composition", "-"*13, "\n")
        print ("{0} = {1} \n {2} = {3} \n {4}={5} \n {6}={7} \n".format(Componentes_f1.value, self.zi[0], Componentes_f2.value, self.zi[1], Componentes_f3.value, self.zi[2], Componentes_f4.value, self.zi[3]))
        print ("Sumatoria zi = {0}".format(np.sum(self.zi)))       
        print ("\n", "-"*13, "Liquid phase composition", "-"*13, "\n")
        print ("{0} = {1} \n {2} = {3} \n {4}={5} \n {6}={7} \n".format(Componentes_f1.value, self.xi[0], Componentes_f2.value, self.xi[1], Componentes_f3.value, self.xi[2], Componentes_f4.value, self.xi[3]))
        print ("Sumatoria xi = {0}".format(np.sum(self.xi)))
        print ("\n", "-"*14, "Vapor phase composition", "-"*13, "\n")
        print ("{0} = {1} \n {2} = {3} \n {4}={5} \n {6}={7} \n".format(Componentes_f1.value, self.yi[0], Componentes_f2.value, self.yi[1], Componentes_f3.value, self.yi[2], Componentes_f4.value, self.yi[3]))
        print ("Sumatoria yi = {0}".format(np.sum(self.yi)))
        print ("-"*53, "\n","Beta = {0}".format(self.Bini), "\n")
        print ("\n","Función R&R = {0}".format(Eg[0]), "\n")
        print ("\n","Derivada función R&R = {0}".format(Eg[1]), "\n", "-"*53)


        return #Eg[0], Eg[1], self.Bini
    
class FlashHP(Fugacidad, Flash):

    def __init__(self, zF):
        Fugacidad.__init__(self, eq, w, Tc, Pc, Tr, R, ep, ni, nT, nC, V, T, P, kij, lij, delta_1, k, Avsl)
        self.zF = zF
        
      
    
    def flash_PT(self):
        # Solución del flash (T,P,ni) isotermico para Ki_(T,P,ni)
        flashID = self.flash_ideal()
        print ("flash (P, T, zi)")
        print ("g, dg, B = ", flashID)
        print ("-"*66)

        self.Bini = flashID[2]
        print ("Beta_r ini = ", self.Bini)
        moles = self.composicion_xy()

        self.xi, self.yi = moles[0], moles[1]
        nil, niv = moles[2], moles[3]

        fi_F = self.fugac()        

        self.Ki = fi_F[0] / fi_F[1]

        L = 1.0

        self.Ki = self.Ki * L

        Ki_1 = self.Ki
        print ("Ki_(P, T, ni) primera = ", self.Ki)

        print ("-"*66)

        #self.Ki = np.array([1.729, 0.832, 0.640])

        #self.Ki = self.wilson(self.Pc, self.Tc, self.w, self.T)
        #print "Ki_(P, T) = ", self.Ki

        while 1:
            i, s = 0, 0.1

            while 1:
                Eg = self.rice()
                print (Eg)
                self.Bini = self.Bini - s * Eg[0] / Eg[1]
                print (self.Bini)
                errorEq = abs(Eg[0])
                i += 1
                #print i

                #if self. Bini < 0 or self.Bini > 1:
                    #break
                #    self.Bini = 0.5
                if i >= 50:
                    pass
                    #break
                if errorEq < 1e-5:
                    break

            print ("Resultado Real = ", Eg)
            print (" Beta r = ", self.Bini)

            moles = self.composicion_xy(zi, self.Ki, self.Bini)
            self.xi, self.yi = moles[0], moles[1]

            #xy = self.composicion_xy(zi, self.Ki, self.Bini)

            print ("C1 -i-C4 n-C4")
            print ("-"*13, "Composición de fase líquida", "-"*13)
            print ("xi = ", moles[0])
            print ("Sxi = ", np.sum(moles[0]))
            print ("-"*13, "Composición de fase vapor", "-"*13)
            print ("yi = ", moles[1])
            print ("Syi = ", np.sum(moles[1]))

            fi_F = self.fugac()

            self.Ki = fi_F[0] / fi_F[1]
            Ki_2 = self.Ki
            dKi = abs(Ki_1 - Ki_2)
            Ki_1 = Ki_2
            print ("Ki_(P, T, ni) = ", self.Ki)

            fun_Ki = np.sum(dKi)
            print ("fun_Ki = ", fun_Ki)

            if fun_Ki < 1e-5:
                break

        return flashID

url = 'Lectura Juan.xlsx'

class DataGPEC():
    
    def __init__(self, url):
        self.url = url
        
    def leerGPEC_1(self):
        """
        El siguiente script python, se puede mejorar generalizando la lectura de etiquetas,
        mientras se pasa la transición GPEC librería
        """
        marcas = ['VAP', 'CRI', 'CEP']
        
        GPEC = pd.read_excel(url)
        
        """
        Revisar las etiquetas, nombre, roturlos de las figurar generadas con este script Python
        para que sean acordes a las variables que se desean gráficar, mientras se automatiza este
        proceso.
        """
        
        
        #------------------------------------------------------------------------------
        DatosGPEC = pd.DataFrame(GPEC)
        VAP = DatosGPEC.loc[(DatosGPEC['T(K)'] == marcas[0])]
        etiquetaVAP = VAP.index.get_values()
        inicioVAP = etiquetaVAP[0]+1
        finalVAP = etiquetaVAP[1]-2
        
        #------------------------------------------------------------------------------
        self.TemperaturaVAP = np.array([DatosGPEC.ix[inicioVAP:finalVAP,0]], dtype=np.float)
        self.PresionVAP = np.array([DatosGPEC.ix[inicioVAP:finalVAP,1]], dtype=np.float)
        self.VolumenLiqVAP = np.array([DatosGPEC.ix[inicioVAP:finalVAP,2]], dtype=np.float)
        self.VolumenVapVAP = np.array([DatosGPEC.ix[inicioVAP:finalVAP,3]], dtype=np.float)
        #------------------------------------------------------------------------------
        CRI = DatosGPEC.loc[(DatosGPEC['T(K)'] == marcas[1])]
        etiquetaCRI = CRI.index.get_values()
        inicioCRI = etiquetaCRI[0]+1
        finalCRI = etiquetaCRI[1]-2
        #------------------------------------------------------------------------------
        self.TemperaturaCRI = np.array([DatosGPEC.ix[inicioCRI:finalCRI,0]], dtype=np.float)
        self.PresionCRI = np.array([DatosGPEC.ix[inicioCRI:finalCRI,1]], dtype=np.float)
        self.VolumenLiqCRI = np.array([DatosGPEC.ix[inicioCRI:finalCRI,2]], dtype=np.float)
        self.VolumenVapCRI = np.array([DatosGPEC.ix[inicioCRI:finalCRI,3]], dtype=np.float)
        #------------------------------------------------------------------------------
        """
        En la segunda línea critica se tiene como referencia el final de la primera línea critica
        y la etiqueta CEP
        """
        
        CEP = DatosGPEC.loc[(DatosGPEC['T(K)'] == marcas[2])]
        etiquetaCEP = CEP.index.get_values()
        inicioCRI_2 = etiquetaCRI[1]+1
        finalCRI_2 = etiquetaCEP[0]-2
        
        self.TemperaturaCRI_2 = np.array([DatosGPEC.ix[inicioCRI_2:finalCRI_2,0]], dtype=np.float)
        self.PresionCRI_2 = np.array([DatosGPEC.ix[inicioCRI_2:finalCRI_2,1]], dtype=np.float)
        self.VolumenLiqCRI_2 = np.array([DatosGPEC.ix[inicioCRI_2:finalCRI_2,2]], dtype=np.float)
        self.VolumenVapCRI_2 = np.array([DatosGPEC.ix[inicioCRI_2:finalCRI_2,3]], dtype=np.float)
        
        return self.TemperaturaCRI_2
    
    def presionVapor(self):
        clear_output()
        pyplot.close("all")
        pyplot.scatter(self.TemperaturaVAP,self.PresionVAP, color = 'red', label = 'Presión de Vapor')
        pyplot.title('Temperatura-Presión')
        pyplot.legend(loc="upper left") 
        pyplot.xlabel('Temperatura [=] K')
        pyplot.ylabel('Presión [=] bar')
        
    def densidadPresion(self):
        clear_output()
        pyplot.close("all")
        pyplot.scatter(self.VolumenLiqVAP,self.PresionVAP, color = 'red', label = 'Líquido')
        pyplot.scatter(self.VolumenVapVAP,self.PresionVAP, color = 'blue', label = 'Vapor')
        pyplot.title('Diagrama Densidad-Presión')
        pyplot.legend(loc="upper right") 
        pyplot.xlabel('Densidad [=] -')  
        pyplot.ylabel('Presión [=] bar')
        
    def diagramaTPcritico(self):
        clear_output()
        pyplot.close("all")
        pyplot.scatter(self.TemperaturaCRI,self.PresionCRI, color = 'red', label = 'Presión Critica')
        pyplot.title('Diagrama Temperatura Cri-Presión Cri')
        pyplot.legend(loc="upper left") 
        pyplot.xlabel('Temperatura [=] K')  
        pyplot.ylabel('Presión [=] bar')
        
    def diagramaDensidadCri(self):
        clear_output()
        pyplot.close("all")
        pyplot.scatter(self.VolumenLiqCRI,self.PresionCRI, color = 'red', label = 'Líquido')
        pyplot.scatter(self.VolumenVapCRI,self.PresionCRI, color = 'blue', label = 'Vapor')
        pyplot.title('Diagrama Densidad Critica')
        pyplot.legend(loc="upper right") 
        pyplot.xlabel('Densidad [=] -')  
        pyplot.ylabel('Presión [=] bar')
        
    def diagramaCritico_2(self):
        clear_output()
        pyplot.close("all")
        fig_2= pyplot.scatter(self.TemperaturaCRI_2,self.PresionCRI_2)
        pyplot.scatter(self.TemperaturaCRI_2,self.PresionCRI_2, color = 'red', label = 'Presión de Critica 2')
        pyplot.title('Diagrama Critico 2')
        pyplot.legend(loc="upper left") 
        pyplot.xlabel('Temperatura [=] K')  
        pyplot.ylabel('Presión [=] bar')
#------------------------------------------------------------------------------

Interfaz "gráfica"


In [12]:
Componentes_1 = widgets.SelectMultiple(
    description="Component 1",
    options=list(Etiquetas))

Componentes_2 = widgets.SelectMultiple(
    description="Component 2",
    options=list(Etiquetas))

button = widgets.Button(description="Upload Data")

def cargarDatos(b):
    clear_output()
    print("Component 1: ", Componentes_1.value)
    Nombre = Componentes_1.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_1 = Propiedades[0]
    Temperatura_Critica_1 = Propiedades[1]
    Presion_Critica_1 = Propiedades[2]
    Z_Critico_1 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor = ", Factor_Acentrico_1)
    print ("Critical Temperature = ", Temperatura_Critica_1, "K")
    print ("Critical Pressure = ", Presion_Critica_1, "bar")
    print ("Z_Critical = ", Z_Critico_1, "\n")
    
    
    print("Component 2: ", Componentes_2.value)
    Nombre = Componentes_2.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_2 = Propiedades[0]
    Temperatura_Critica_2 = Propiedades[1]
    Presion_Critica_2 = Propiedades[2]
    Z_Critico_2 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor  = ", Factor_Acentrico_2)
    print ("Critical Temperature = ", Temperatura_Critica_2, "K")
    print ("Critical Pressure = ", Presion_Critica_2, "bar")
    print ("Z_Critical = ", Z_Critico_2)
    
    global TcDato, PcDato, wDato
    
    TcDato = np.array([Temperatura_Critica_1, Temperatura_Critica_2])
    PcDato = np.array([Presion_Critica_1, Presion_Critica_2])
    wDato = np.array([Factor_Acentrico_1, Factor_Acentrico_2])


button.on_click(cargarDatos)
#display(button)

page1 = widgets.VBox(children=[Componentes_1, Componentes_2, button], padding=4)


#VBox([VBox([Button(description='Press'), Dropdown(options=['a', 'b']), Button(description='Button')]), 
#      VBox([Button(), Checkbox(), IntText()]), 
#      VBox([Button(), IntSlider(), Button()])], background_color='#EEE')


ecuacionEstado = widgets.Dropdown(description='Fluid :', padding=4, options=['SRK', 'PR', 'RKPR'])
modeloSolido = widgets.Dropdown(description='Solid :', padding=4, options=['Model I', 'Model II', 'Model III'])

button = widgets.Button(description="Upload Models")

def cargarModelos(b):
    clear_output()
    global eq    
    eq = ecuacionEstado.value
    
    print("Component 1: ", Componentes_1.value)
    print("Component 2: ", Componentes_2.value)
    
    print("Fluid Model : ", ecuacionEstado.value)
    print("Solid Model : ", modeloSolido.value)
    
    


button.on_click(cargarModelos)

page2 = widgets.Box(children=[ecuacionEstado, modeloSolido, button], padding=4)


Temp_ini = widgets.Text(description='Initial', padding=4, value="0.0")
Temp_fin = widgets.Text(description='Final', padding=4, value="0.0")

Pres_ini = widgets.Text(description='Initial', padding=4, value="0.0")
Pres_fin = widgets.Text(description='Final', padding=4, value="0.0")

n1 = widgets.Text(description='Mole light component', padding=4, value="0.0")
n2 = widgets.Text(description='Mole heavy component', padding=4, value="0.0")

#button = widgets.Button(description="Cargar Condiciones")

titulo = widgets.HTML(value="<C><H1> System Conditions <H1>")
tempe_info = widgets.HTML(value="<C><H3> Temperature <H3>")
press_info = widgets.HTML(value="<C><H3> Pressure <H3>")
fluid_info = widgets.HTML(value="<C><H3> Mole fracction in the fluid <H3>")

button = widgets.Button(description="Upload Conditions")

def cargarParametros(b):
    clear_output()
    
    global initial_temperature, initial_pressure, nif
    
    initial_temperature = float(Temp_ini.value)
    initial_pressure = float(Pres_ini.value)
    nif = np.array([float(n1.value), float(n2.value)])
        
    print("Component 1: ", Componentes_1.value)
    print("Component 2: ", Componentes_2.value)
    
    print("Fluid Model : ", ecuacionEstado.value)
    print("solid Model : ", modeloSolido.value)
    
    print("Initial_temperature = ", initial_temperature, type(initial_temperature))
    print("Final_temperature = ", Temp_fin.value)
    
    print("Initial_pressure =", initial_pressure, type(initial_pressure))
    print("Final_pressure =", Pres_fin.value)
    
    print("Mole fraccion light component n1 =", n1.value)
    print("Mole fraccion heavy component n2 =", n2.value)
    
    print("Mole fracction in the fluid ", nif) 
    
    print(initial_temperature, type(initial_temperature))
    

button.on_click(cargarParametros)

page3 = widgets.Box(children=[titulo, tempe_info, Temp_ini, Temp_fin, press_info, Pres_ini, Pres_fin, fluid_info, n1, n2, button], padding=4)






button = widgets.Button(description="Solid-Fluid")
#display(button)


nnCC_1 = 1
nnCC_2 = 2

def calcularSolidoFluido(b):
    clear_output()
    #Tcal = fsolve(equilibrioSF,guess,args=(Pe, nnCC_1, nnCC_2), xtol=1e-4)
    
    #initial_temperature = [346.5] # T [=] K
    #initial_pressure = 137.9 # [=] bar
    
    Tcal = fsolve(equilibrioSF,initial_temperature,args=(initial_pressure, nif, 1, 2, Avsl), xtol=1e-4)
    print("Temperature ESF = ", Tcal, "K")
    

button.on_click(calcularSolidoFluido)
#display(button)


page4 = widgets.Box(children=[button], padding=4)

button = widgets.Button(description="Diagram Solid-Fluid")

def DiagramaSolidoFluido(b):
    clear_output()
    #Tcal = fsolve(equilibrioSF,guess,args=(Pe, 1, 2), xtol=1e-4)
    #Tcal = fsolve(equilibrioSF,guess,args=(Pe, nnCC_1, nnCC_2), xtol=1e-4)
    initial_temperature =346.5 # T [=] K
    initial_pressure = 136.9 # [=] bar
    
    #346.5 136.9
    # n1, n2 = 1, 2 por defecto para el equilibrio sólido-fluido
    Tcal = fsolve(equilibrioSF,initial_temperature,args=(initial_pressure, nif, 1, 2, Avsl), xtol=1e-4)

    print(Tcal, "K")
    
    pyplot.scatter(Tres,pres, color = 'red', label = 'PR')
    pyplot.scatter(temp,pres, label = 'Data')
    pyplot.title('Temperature Equilibrium Solid Liquid')
    pyplot.legend(loc="upper left") 
    pyplot.xlabel('Temperature [=] K')  
    pyplot.ylabel('Pressure [=] bar')
    

button.on_click(DiagramaSolidoFluido)


page5 = widgets.Box(children=[button], padding=4)




DatosTemperatura_Exp = np.array([323.65, 326.04, 326.43])
#DatosTemperatura_Exp = np.array([323.65, 326.04, 326.43, 328.12])

DatosPresionp_Exp = np.array([1.0, 101.0, 136.9])
#DatosPresionp_Exp = np.array([1.0, 101.0, 136.9, 183.8])


posicion = np.arange(len(DatosPresionp_Exp))
TemperaturasModelo = np.ones((len(DatosPresionp_Exp)))
TemperaturasModelo

Avsl = -0.32595074


button = widgets.Button(description="Regression of Parameters")

def regresionParametros(b):
    clear_output()
    
    def minimizarVSL(Avsl):
        for T, P, i in zip(DatosTemperatura_Exp, DatosPresionp_Exp, posicion):
            print ("Initial Temperature = ", T, "K", "Pressure = ", P, "bar", "Experimental Data = ", i+1)
            initial_temperature = T # T [=] K
            initial_pressure = P # [=] bar
            # tol=
            TemperaturasModelo[i] = fsolve(equilibrioSF,initial_temperature,args=(initial_pressure, nif, 1, 2, Avsl), xtol=1e-4)
            
        funcionObjetivo = np.sum((DatosTemperatura_Exp - TemperaturasModelo) ** 2)
        print("modelTemperature = ", TemperaturasModelo)
        print("Objective Function = ", funcionObjetivo)
        
        return funcionObjetivo
    
    opt = sp.optimize.minimize(minimizarVSL, Avsl, method='L-BFGS-B')
    
    print("optimal parameter", opt)


button.on_click(regresionParametros)


page6 = widgets.Box(children=[button], padding=4)

tabs = widgets.Tab(children=[page1, page2, page3, page4, page5, page6])
#display(tabs)

tabs.set_title(0, 'Components')
tabs.set_title(1, 'Models')
tabs.set_title(2, 'Conditions')
tabs.set_title(3, 'Results')
tabs.set_title(4, 'Experimental Data')
tabs.set_title(5, 'Regression of Parameters')

#--------------------- flash Isothermal------------------------------

Componentes_f1 = widgets.SelectMultiple(
    description="Name 1",
    options=list(Etiquetas))

Componentes_f2 = widgets.SelectMultiple(
    description="Name 2",
    options=list(Etiquetas))

Componentes_f3 = widgets.SelectMultiple(
    description="Name 3",
    options=list(Etiquetas))

Componentes_f4 = widgets.SelectMultiple(
    description="Name 4",
    options=list(Etiquetas))

button = widgets.Button(description="Upload Data")

def cargarDatos(b):
    clear_output()
    print("Component 1: ", Componentes_f1.value)
    Nombre = Componentes_f1.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_1 = Propiedades[0]
    Temperatura_Critica_1 = Propiedades[1]
    Presion_Critica_1 = Propiedades[2]
    Z_Critico_1 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor = ", Factor_Acentrico_1)
    print ("Critical Temperature = ", Temperatura_Critica_1, "K")
    print ("Critical Pressure = ", Presion_Critica_1, "bar")
    print ("Z_Critical = ", Z_Critico_1, "\n")
    
    
    print("Component 2: ", Componentes_f2.value)
    Nombre = Componentes_f2.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_2 = Propiedades[0]
    Temperatura_Critica_2 = Propiedades[1]
    Presion_Critica_2 = Propiedades[2]
    Z_Critico_2 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor  = ", Factor_Acentrico_2)
    print ("Critical Temperature = ", Temperatura_Critica_2, "K")
    print ("Critical Pressure = ", Presion_Critica_2, "bar")
    print ("Z_Critical = ", Z_Critico_2, "\n")
    
    print("Component 3: ", Componentes_f3.value)
    Nombre = Componentes_f3.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_3 = Propiedades[0]
    Temperatura_Critica_3 = Propiedades[1]
    Presion_Critica_3 = Propiedades[2]
    Z_Critico_3 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor = ", Factor_Acentrico_3)
    print ("Critical Temperature = ", Temperatura_Critica_3, "K")
    print ("Critical Pressure = ", Presion_Critica_3, "bar")
    print ("Z_Critical = ", Z_Critico_3, "\n")
    
    
    print("Component 4: ", Componentes_f4.value)
    Nombre = Componentes_f4.value
    Propiedades = data2.loc[Nombre]
    Factor_Acentrico_4 = Propiedades[0]
    Temperatura_Critica_4 = Propiedades[1]
    Presion_Critica_4 = Propiedades[2]
    Z_Critico_4 = Propiedades[3]

    #print(Propiedades)
    print ("Acentric Factor  = ", Factor_Acentrico_4)
    print ("Critical Temperature = ", Temperatura_Critica_4, "K")
    print ("Critical Pressure = ", Presion_Critica_4, "bar")
    print ("Z_Critical = ", Z_Critico_4, "\n")
    
    global TcDato_f, PcDato_f, wDato_f
    
    TcDato_f = np.array([Temperatura_Critica_1, Temperatura_Critica_2, Temperatura_Critica_3, Temperatura_Critica_4])
    PcDato_f = np.array([Presion_Critica_1, Presion_Critica_2, Presion_Critica_3, Presion_Critica_4])
    wDato_f = np.array([Factor_Acentrico_1, Factor_Acentrico_2, Factor_Acentrico_3, Factor_Acentrico_4])


button.on_click(cargarDatos)
#display(button)

page_f1 = widgets.VBox(children=[Componentes_f1, Componentes_f2, Componentes_f3, Componentes_f4, button], padding=4)


#------------------ page_f2
ecuacionEstado_f = widgets.Dropdown(description='Fluid :', padding=4, options=['SRK', 'PR', 'RKPR'])

button = widgets.Button(description="Upload Models")

def cargarModelos(b):
    clear_output()
    global eq    
    eq = ecuacionEstado.value
    
    print("Component 1: ", Componentes_f1.value)
    print("Component 2: ", Componentes_f2.value)
    print("Component 3: ", Componentes_f3.value)
    print("Component 4: ", Componentes_f4.value)    
    print("Fluid Model : ", ecuacionEstado_f.value)
    
    
button.on_click(cargarModelos)

page_f2 = widgets.Box(children=[ecuacionEstado_f, button], padding=4)

#------------------ page_f2

#------------------ page_f3
Temp_ini_f = widgets.Text(description='K', padding=4, value="0.0")

Pres_ini_f = widgets.Text(description='Bar', padding=4, value="0.0")

n1_f = widgets.Text(description='Name 1', padding=4, value="0.0")
n2_f = widgets.Text(description='Name 2', padding=4, value="0.0")
n3_f = widgets.Text(description='Name 3', padding=4, value="0.0")
n4_f = widgets.Text(description='Name 4', padding=4, value="0.0")


titulo = widgets.HTML(value="<C><H1> System Conditions <H1>")
tempe_info = widgets.HTML(value="<C><H3> Temperature <H3>")
press_info = widgets.HTML(value="<C><H3> Pressure <H3>")
fluid_info = widgets.HTML(value="<C><H3> Mole fracction in the fluid <H3>")

button = widgets.Button(description="Upload Conditions")

def cargarParametros(b):
    clear_output()
    
    global zi_F, temperature_f, pressure_f, nif
    
    temperature_f = float(Temp_ini_f.value)
    pressure_f = float(Pres_ini_f.value)
    zi_F = np.array([float(n1_f.value), float(n2_f.value), float(n3_f.value), float(n4_f.value)])
    nif = np.array([float(n1_f.value), float(n2_f.value), float(n3_f.value), float(n4_f.value)])
        
    print("Component 1: ", Componentes_f1.value)
    print("Component 2: ", Componentes_f2.value)
    print("Component 3: ", Componentes_f3.value)
    print("Component 4: ", Componentes_f4.value)
    
    print("Fluid Model : ", ecuacionEstado_f.value)
    
    print("Temperature_f = ", temperature_f, type(temperature_f))
    
    print("Pressure_f = ", pressure_f, type(pressure_f))
    
    print("Mole fraccion component 1 = ", n1_f.value)
    print("Mole fraccion component 2 = ", n2_f.value)
    print("Mole fraccion component 3 = ", n3_f.value)
    print("Mole fraccion component 4 = ", n4_f.value)
    
    print("Mole fracction in the fluid = ", zi_F, type(zi_F)) 
    
    print(temperature_f, type(temperature_f))
    

button.on_click(cargarParametros)

page_f3 = widgets.Box(children=[titulo, tempe_info, Temp_ini_f, press_info, Pres_ini_f, fluid_info, n1_f, n2_f, n3_f, n4_f, button], padding=4)
#------------------ page_f3

#------------------ page_f4

button = widgets.Button(description="Flash Calculation")

def calcularFlashPT(b):
    clear_output()
    #Tcal = fsolve(equilibrioSF,guess,args=(Pe, nnCC_1, nnCC_2), xtol=1e-4)
    
    #initial_temperature = [346.5] # T [=] K
    #initial_pressure = 137.9 # [=] bar
    fhid = Flash(zi_F, temperature_f, pressure_f, TcDato_f, PcDato_f, wDato_f)
    fhid.flash_ideal()
    

button.on_click(calcularFlashPT)
#display(button)


page_f4 = widgets.Box(children=[button], padding=4)



#------------------ page_f4
flash = widgets.Tab(children=[page_f1, page_f2, page_f3, page_f4])

flash.set_title(0, 'Components')
flash.set_title(1, 'Models')
flash.set_title(2, 'Conditions')
flash.set_title(3, 'Results')
#tabs.set_title(4, 'Experimental Data')
#tabs.set_title(5, 'Regression of Parameters')

#--------------------- GPEC ------------------------------

name_GPEC = widgets.Text(description='File name', padding=4, value=" ")
url = name_GPEC.value

titulo = widgets.HTML(value="<C><H1> Data GPEC <H1>")

button_1 = widgets.Button(description="UpData GPEC")

def upGPEC(b):
    clear_output()
    
    DGPEC = DataGPEC(url)
    DGPEC.leerGPEC_1()
    print ("Upload {0}".format(url))
        

button_1.on_click(upGPEC)

button_2 = widgets.Button(description="Vapor pressure")

def diagram_1(b):
    clear_output()
    DGPEC = DataGPEC(url)
    DGPEC.leerGPEC_1()
    DGPEC.presionVapor()
    
button_2.on_click(diagram_1)


button_3 = widgets.Button(description="Diagram Density-Pressure")

def diagram_2(b):
    clear_output()
    DGPEC = DataGPEC(url)
    DGPEC.leerGPEC_1()
    DGPEC.densidadPresion()
    

button_3.on_click(diagram_2)

page_G1 = widgets.Box(children=[titulo, name_GPEC, button_1, button_2, button_3], padding=4)

gpec = widgets.Tab(children=[page_G1])

gpec.set_title(0, 'Upload Data')



accord = widgets.Accordion(children=[tabs, flash, gpec], width=400)
display(accord)




accord.set_title(0, 'Pure Solid-Binary Fluid')
accord.set_title(1, 'Isothermal Flash Calculation')
accord.set_title(2, 'Data GPEC')
#accord.set_title(3, 'Regression of Parameters Solid-Fluid')
#accord.set_title(4, 'Pure Fluid')



# 346.5 136.9

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [13]:
url = 'Lectura Juan.xlsx'