Flash_calculation.py

import math
import numpy as np
from scipy import optimize
import matplotlib.pyplot as plt
import mpltern
from PT_flash import PTFlash
from PH_flash import PHFlash

class FlashCalculation:
    """
     This part is dedicated for running all routines attached to Flash calculations.
     In this class the following methods are present:
        - flash_calculation: Calculation of flash for all compositions and print the results
        - f_obj: Cost function for adjustment of Kij, in case of unknown such parameter previously
        - Kij_adjustment: Optimization of flash function for obtaining Kij parameter
        - graph: Generation of P-xy or T-xy plot for binary substances
        - ARD: Average Relative Deviation between experimental and calcutated compositions
        - ternary_plot: Generation of ternary diagram
    """

    def __init__(self):
        pass

    @staticmethod
    def flash_calculation(P, R, T, z, Tc, Pc, w, Kij, P_sat, Hf, T_range, A, B, C, D, sigma, Ek, mi, Components, type, Eos):

        equations = ['SRK', 'Peng-Robinson', 'PC-SAFT']

        if Eos >= len(equations):
            raise TypeError("This equation was not implemented. Please check the available EoS:\n"
                            "- SRK ---------------> EoS=0\n"
                            "- Peng-Robinson -----> EoS=1\n"
                            "- PC-SAFT -----------> EoS=2")
        else:
            if type == 1:
                Ki = (Pc/P)*np.exp(5.37*(1+w)*(1-(Tc/T)))
                F, V, x, y, n_ite, Ki_new = PTFlash().flash(P, R, T, z, Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                                mi, A, B, C, D, equation=equations[Eos])
                print("")
                print(f'################ PT flash accomplished successfully ##################')
                print(f'Selected equation of state: {equations[Eos]}')
                print(f"The number of iterations: {n_ite}")
                print(f"Volume fraction in the liquid phase: {F}")
                print(f"Volume fraction in the vapor phase {V}")
                for i, comp in enumerate(Components):
                    print(f'Molar fraction in the liquid phase for {comp} is {x[i]:.3f} in equilibrium')
                    print(f'Molar fraction in the vapor phase for {comp} is {y[i]:.3f} in equilibrium')

            elif type == 2:
                n_ite, T, x, y, V, F = PHFlash().flash(Hf, P[0], R, T_range, Tc, z, Pc, w, Kij, P_sat, sigma, Ek,
                                                   mi, A, B, C, D, equation=equations[Eos])
                print("")
                print(f'################ PH flash accomplished successfully ##################')
                print(f'Selected equation of state: {equations[Eos]}')
                print(f"The number of iterations: {n_ite}")
                print(f'Tank temperature {T:.2f} K')
                print(f"Volume fraction in the liquid phase: {F:.2f} ")
                print(f"Volume fraction in the vapor phase {V:.2f}")
                for i, comp in enumerate(Components):
                    print(f'Molar fraction in the liquid phase for {comp} is {x[i]:.2f} in equilibrium')
                    print(f'Molar fraction in the vapor phase for {comp} is {y[i]:.2f} in equilibrium')

            else:
                raise TypeError(f"You selected type {type}. Such option doesn't exist. The available options are:\n"
                                "- PT flash --------> type=1\n"
                                "- PH flash --------> type=2")

    @staticmethod
    def f_obj(Kij, X_exp, Y_exp, z, P, R, T, Tc, Pc, w, P_sat, A, B, C, D, sigma, Ek, mi, Eos):

        cont = 1
        equations = ['SRK', 'Peng-Robinson', 'PC-SAFT']

        #For P-x-y
        if len(P)>1:
            # Initial guess of Ki (Wilson equation)
            Ki     = (Pc/P[0])*np.exp(5.37*(1+w)*(1-(Tc/T)))
            X      = np.zeros(len(P))
            Y      = np.zeros(len(P))
            x0_exp = np.zeros(len(P))
            y0_exp = np.zeros(len(P))
            for i in range(len(P)):
                _, _, x, y, _, Ki_new = PTFlash().flash(P[i], R, T, z[i], Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                                    mi, A, B, C, D, equation=equations[Eos])
                Ki        = Ki_new
                X[i]      = x[0]
                Y[i]      = y[0]
                x0_exp[i] = X_exp[i][0]
                y0_exp[i] = Y_exp[i][0]

        #For T-x-y
        elif len(T)>1:
            # Initial guess of Ki (Wilson equation)
            Ki     = (Pc/P)*np.exp(5.37*(1+w)*(1-(Tc/T[0])))
            X      = np.zeros(len(T))
            Y      = np.zeros(len(T))
            x0_exp = np.zeros(len(T))
            y0_exp = np.zeros(len(T))
            for i in range(len(T)):
                _, _, x, y, _, Ki_new = PTFlash().flash(P[0], R, T[i], z[i], Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                                mi, A, B, C, D, equation=equations[Eos])
                Ki = Ki_new
                X[i] = x[0]
                Y[i] = y[0]
                x0_exp[i] = X_exp[i][0]
                y0_exp[i] = Y_exp[i][0]

        Pa = abs(X-x0_exp)
        Pb = abs(Y-y0_exp)

        f_obj = (1/len(x0_exp))*sum(Pa+Pb)

        return f_obj

    @classmethod
    def Kij_adjustment(cls, X_exp, Y_exp, z, P, R, T, Tc, Pc, w, P_sat, A, B, C, D, sigma, Ek, mi, Eos):
        """
        * Binary interaction parameter estimation by Nelder-Mead optimization algorithm
        """

        args   = (X_exp, Y_exp, z, P, R, T, Tc, Pc, w, P_sat, A, B, C, D, sigma, Ek, mi, Eos)
        n      = math.factorial(len(w))/(math.factorial(2)*(math.factorial(len(w)-2)))
        x0     = np.zeros(int(n))
        x0 = np.array([0.09046009537009686])
        bounds = []
        bound  = (0.0, 0.9)
        for i in range(len(x0)):
            bounds.append(bound)

        res = optimize.minimize(cls.f_obj, x0, args, method='Nelder-Mead', options={'fatol':1e-10}, bounds=tuple(bounds))
        print("")
        print(f'################ Kij adjustment accomplished successfully ##################')
        for i, kij in enumerate(res.x):
            print(f'Kij{i}: {kij:.4f}')

        return res.x

    @classmethod
    def binary_plot(cls, Components, x_exp, y_exp, z, Pressure, R, Temperature, Tc, Pc, w, Kij, P_sat, A, B, C, D, sigma, Ek, mi, Eos, use_Kij, save=0):
        """
        * Plot of adjusted model
        """
        # Other general parameters
        equations = ['SRK', 'Peng-Robinson', 'PC-SAFT']

        colors  = ['r','b','k','c','m','y','g','gray','purple','cyan','orange']

        # Plot general parameters
        fig, ax = plt.subplots(1, 1, figsize=(10.8, 8.8))

        # ----------------------- Pxy DIAGRAM ---------------------- #
        if type(Pressure) is dict:
            NAME = f"{Components[0] + '_' + Components[1] + '_' + equations[Eos]+ '_' + str(Temperature)}"
            for k, T in enumerate(Temperature):
                c = colors[k]

                if Eos == 0:
                    legend_model = f'{"SRK "+str(T)+"K"}'
                elif Eos == 1:
                    legend_model = f'{"PR "+ str(T)+"K"}'
                elif Eos == 2:
                    legend_model = f'{"PC-SAFT "+str(T)+"K"}'

                if use_Kij == 0:
                    Kij = cls.Kij_adjustment(x_exp[T], y_exp[T], z[T], Pressure[T], R, T, Tc, Pc, w, P_sat, A, B, C, D, sigma, Ek, mi, Eos)

                # Initial guess of Ki (Wilson equation)
                Ki    = (Pc/Pressure[T][0])*np.exp(5.37*(1+w)*(1-(Tc/T)))
                x_out = np.zeros(len(Pressure[T]))
                y_out = np.zeros(len(Pressure[T]))
                X_exp, Y_exp = [], []
                y_label = "Pressure [MPa]"

                for i in range(len(Pressure[T])):
                    _, _, x, y, _, Ki_new = PTFlash().flash(Pressure[T][i], R, T, z[T][i], Tc, Pc, w, Ki, Kij[T], P_sat, sigma, Ek,
                                                    mi, A, B, C, D, equation=equations[Eos])
                    Ki = Ki_new
                    x_out[i] = x[0]
                    y_out[i] = y[0]
                    X_exp.append(x_exp[T][i][0])
                    Y_exp.append(y_exp[T][i][0])

                if use_Kij == 0:
                    print("")
                    print(f'################ Kij adjustment accomplished successfully ##################')
                    for j, kij in enumerate(Kij):
                        print(f'Kij{j}: {kij}')


                if k == 0:
                    ax.plot(X_exp, Pressure[T]/1e6, 'ko', label='Dados experimentais', markersize=10)
                else:
                    ax.plot(X_exp, Pressure[T]/1e6, 'ko', markersize=10)

                ax.plot(Y_exp, Pressure[T]/1e6, 'ko', markersize=10)
                ax.plot(x_out, Pressure[T]/1e6, c, label=legend_model, linewidth=2.5)
                ax.plot(y_out, Pressure[T]/1e6, c, linewidth=2.5)
                #ax.set_xlabel(f"x\u2081,y\u2081", fontsize=25)
                ax.set_xlabel(f"x\u2081,y\u2081", fontsize=30)
                ax.set_ylabel(y_label, fontsize=30)
                ax.tick_params(labelsize=30)
                plt.gca().spines['bottom'].set_linewidth(4)
                plt.gca().spines['left'].set_linewidth(4)
                plt.gca().spines['top'].set_linewidth(4)
                plt.gca().spines['right'].set_linewidth(4)
                plt.legend(loc='upper left', fontsize=22, frameon=False)
            if save == 1:
                plt.savefig(f"Output/{NAME}.png", dpi=300)

        # ----------------------- Txy DIAGRAM ---------------------- #
        elif type(Temperature) is dict:
            NAME = f"{Components[0] + '_' + Components[1] + '_' + equations[Eos] + '_' + str(Pressure/1e6) + 'MPa'}"
            for i, P in enumerate(Pressure):

                if Eos == 0:
                    legend_model = f'{"SRK EdE" + "_" + str(P)}'
                elif Eos == 1:
                    legend_model = f'{"PR EdE" + "_" + str(P)}'
                elif Eos == 2:
                    legend_model = f'{"PC-SAFT EdE" + "_" + str(P)}'

                if use_Kij == 0:
                    Kij = cls.Kij_adjustment(x_exp[P], y_exp[P], z[P], Pressure, R, Temperature[P], Tc, Pc, w, P_sat, A, B, C, D, sigma, Ek, mi, Eos)

                # Initial guess of Ki (Wilson equation)
                Ki = (Pc/P)*np.exp(5.37*(1+w)*(1-(Tc/Temperature[P][0])))

                x_out = np.zeros(len(Temperature[P]))
                y_out = np.zeros(len(Temperature[P]))
                X_exp, Y_exp = [], []
                y_label = "Temperature (K)"
                for i in range(len(Temperature[P])):

                    _, _, x, y, _, Ki_new = PTFlash().flash(P, R, Temperature[P][i], z[P][i], Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                                    mi, A, B, C, D, equation=equations[Eos])
                    Ki = Ki_new
                    x_out[i] = x[0]
                    y_out[i] = y[0]
                    X_exp.append(x_exp[P][i][0])
                    Y_exp.append(y_exp[P][i][0])

                if use_Kij == 0:
                    print("")
                    print(f'################ Kij adjustment accomplished successfully ##################')
                    for i, kij in enumerate(Kij):
                        print(f'Kij{i}: {kij}')

                ax.plot(X_exp, Temperature[P], 'ko', label='Experimental data', markersize=10)
                ax.plot(Y_exp, Temperature[P], 'ko', markersize=10)
                ax.plot(x_out, Temperature[P], 'black', label=legend_model, linewidth=2.5)
                ax.plot(y_out, Temperature[P], 'black', linewidth=2.5)
                ax.set_xlabel(f"x\u2081,y\u2081", fontsize=30)
                ax.set_ylabel(y_label, fontsize=30)
                ax.tick_params(labelsize=30)
                plt.gca().spines['bottom'].set_linewidth(4)
                plt.gca().spines['left'].set_linewidth(4)
                plt.gca().spines['top'].set_linewidth(4)
                plt.gca().spines['right'].set_linewidth(4)
                plt.legend(loc='best', fontsize=16, frameon=False)

            if save == 1:
                plt.savefig(f"Output/{NAME}.png", dpi=300)

        plt.show()

    @classmethod
    def AARD(cls, Components, X_exp, Y_exp, Z_exp, P, R, T, Tc, Pc, w, Kij, P_sat, A, B, C, D, sigma, Ek, mi, Eos, save):

        equations = ['SRK', 'Peng-Robinson', 'PC-SAFT']
        X     = np.zeros(len(X_exp[T]))
        Y     = np.zeros(len(Y_exp[T]))
        x_exp = np.zeros(len(X_exp[T]))
        y_exp = np.zeros(len(Y_exp[T]))
        z     = Z_exp

        # For P-x-y
        if type(P) is dict:
            # Initial guess of Ki (Wilson equation)
            Ki = (Pc/P[T][0])*np.exp(5.37* (1+w)*(1-(Tc/T)))
            for i in range(len(P[T])):
                _, _, x, y, _, Ki_new = PTFlash().flash(P[T][i], R, T, z[T][i], Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                            mi, A, B, C, D, equation=equations[Eos])
                Ki   = Ki_new
                X[i] = x[0]
                Y[i] = y[0]
                x_exp[i] = X_exp[T][i][0]
                y_exp[i] = Y_exp[T][i][0]

        # For T-x-y
        elif type(T) is dict:
            # Initial guess of Ki (Wilson equation)
            Ki = (Pc/P)*np.exp(5.37*(1+w)*(1-(Tc/T[P][0])))
            for i in range(len(T[P])):
                _, _, x, y, _, Ki_new = PTFlash().flash(P, R, T[P][i], z[P][i], Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                                        mi, A, B, C, D, equation=equations[Eos])
                Ki   = Ki_new
                X[i] = x[0]
                Y[i] = y[0]
                x_exp[i] = X_exp[T][i][0]
                y_exp[i] = Y_exp[T][i][0]

        x = X[1:]
        y = Y[1:]
        ARD = (100 / len(x)) * (sum(abs((x - x_exp[1:]) / x_exp[1:])) + sum(abs((y - y_exp[1:]) / y_exp[1:])))

        print(f'\n ######################### FINISHED ############################\n'
              f'ARD: {ARD:.2f}')

        if save == 1:
            import pandas as pd
            NAME = f"{Components[0]+'_'+Components[1]+'_'+equations[Eos]+'_'+str(T[0])+'_'+'K'}"
            ARD_list = np.array(len(X_exp)*[ARD])
            T_list   = np.array(len(X_exp)*[T[0]])
            output = {'x': X,
                    'y': Y,
                    'x_exp': X_exp,
                    'y_exp': Y_exp,
                    'P': P,
                    'T': T_list,
                    'ARD': ARD_list}

            df = pd.DataFrame(output)
            df.to_csv(f'Output/{NAME}.csv', index=False)

    @classmethod
    def ternary_plot(cls, P, R, T, z, x_exp, y_exp, Tc, Pc, w, Kij, P_sat, A, B, C, D, sigma, Ek, mi, Components, Eos, save):

        equations = ['SRK', 'Peng-Robinson', 'PC-SAFT']

        if Eos >= len(equations):
            raise TypeError("This equation was not implemented. Please check the available EoS:\n"
                            "- SRK ---------------> EoS=0\n"
                            "- Peng-Robinson -----> EoS=1\n"
                            "- PC-SAFT -----------> EoS=2")
        else:


            if len(P) != 0:
                data_storage_x_EoS = []
                data_storage_y_EoS = []
                data_storage_x_exp = []
                data_storage_y_exp = []

                for j, pressure in enumerate(P):
                    # Initial guess of Ki (Wilson equation)
                    Ki = (Pc/P[0])*np.exp(5.37*(1+w)*(1-(Tc/T)))
                    x1_output, x2_output, x3_output = [], [], []
                    y1_output, y2_output, y3_output = [], [], []
                    x1_exp, x2_exp, x3_exp = [], [], []
                    y1_exp, y2_exp, y3_exp = [], [], []
                    for i in range(len(z[pressure])):
                        F, V, x, y, n_ite, Ki_new = PTFlash().flash(P[j], R, T[0], z[pressure][i], Tc, Pc, w, Ki, Kij, P_sat, sigma, Ek,
                                                            mi, A, B, C, D, equation=equations[Eos])
                        Ki = Ki_new
                        x1_output.append(x[0])
                        x2_output.append(x[1])
                        x3_output.append(x[2])

                        y1_output.append(y[0])
                        y2_output.append(y[1])
                        y3_output.append(y[2])

                        # Organize the experimental points
                    for i in range(len(x_exp[pressure])):
                        x1_exp.append(x_exp[pressure][i][0])
                        x2_exp.append(x_exp[pressure][i][1])
                        x3_exp.append(x_exp[pressure][i][2])

                        y1_exp.append(y_exp[pressure][i][0])
                        y2_exp.append(y_exp[pressure][i][1])
                        y3_exp.append(y_exp[pressure][i][2])

                    data_x_EoS = {Components[0]: x1_output, Components[1]: x2_output, Components[2]: x3_output}
                    data_x_exp = {Components[0]: x1_exp, Components[1]: x2_exp, Components[2]: x3_exp}
                    data_y_EoS = {Components[0]: y1_output, Components[1]: y2_output, Components[2]: y3_output}
                    data_y_exp = {Components[0]: y1_exp, Components[1]: y2_exp, Components[2]: y3_exp}

                    data_storage_x_EoS.append(data_x_EoS)
                    data_storage_y_EoS.append(data_y_EoS)
                    data_storage_x_exp.append(data_x_exp)
                    data_storage_y_exp.append(data_y_exp)

                # ----------------------------- Ternary diagram ---------------------------#
                fig = plt.figure(figsize=(10.8, 8.8))
                ax  = fig.add_subplot(projection="ternary")
                ax.set_tlabel(Components[1], fontsize=18)
                ax.set_llabel(Components[0], fontsize=18)
                ax.set_rlabel(Components[2], fontsize=18)
                ax.tick_params(labelsize=16)

                markers = (len(P)+2)*['o', '^']
                colors  = ['k','r','b','c','m','y','g','gray','purple','cyan','orange']
                n = 0
                for i in range(len(data_storage_x_EoS)):
                    data_x_EoS = data_storage_x_EoS[i]
                    data_x_exp = data_storage_x_exp[i]
                    data_y_EoS = data_storage_y_EoS[i]
                    data_y_exp = data_storage_y_exp[i]
                    ax.plot(data_x_EoS[Components[1]], data_x_EoS[Components[0]], data_x_EoS[Components[2]], colors[i], label=f'{round(P[i]/10e5,2)} MPa')
                    ax.plot(data_x_exp[Components[1]], data_x_exp[Components[0]], data_x_exp[Components[2]], colors[i]+markers[i+n])
                    ax.plot(data_y_EoS[Components[1]], data_y_EoS[Components[0]], data_y_EoS[Components[2]], colors[i])
                    ax.plot(data_y_exp[Components[1]], data_y_exp[Components[0]], data_y_exp[Components[2]], colors[i]+markers[i+n+1])
                    n += 1

                ax.grid()
                a, = ax.plot(-1.9, 1.9, 1.9, '^k', markerfacecolor='white')
                b, = ax.plot(-1.9, 1.9, 1.9, 'ok', markerfacecolor='white')
                legend1 = plt.legend([a, b], ['Vapor phase', 'Liquid phase'], title='Experimental data', loc='upper left', edgecolor='0.0', title_fontsize=16, fontsize=14)
                plt.gca().add_artist(legend1)
                plt.legend(title='Pressure', loc='best', edgecolor='0.0', title_fontsize=16, fontsize=14)
                NAME = f"{Components[0] + '_' + Components[1] + '_' + Components[1] + '_' + equations[Eos] + '_' + str(T[0]) + '_' + 'K'}"

                if save == 1:
                    plt.savefig(f'Output/{NAME}.png', dpi=300)

                plt.show()

    @classmethod
    def graph(cls, P, R, T, z, x_exp, y_exp, Tc, Pc, w, Kij, P_sat, A, B, C, D, sigma, Ek, mi, Components, Eos, use_Kij, save=0):

        keys = list(x_exp.keys())

        if len(x_exp[keys[0]][0]) == 2:
            cls.binary_plot(Components, x_exp, y_exp, z, P, R, T, Tc, Pc, w, Kij, P_sat,
                            A, B, C, D, sigma, Ek, mi, Eos, use_Kij, save)

        elif len(x_exp[keys[0]][0]) == 3:
            cls.ternary_plot(P, R, T, z, x_exp, y_exp, Tc, Pc, w, Kij, P_sat,
                             A, B, C, D, sigma, Ek, mi, Components, Eos, save)

        else:
            raise TypeError("Check how many components exist. The available number of components for "
                            "generating the diagrams are 2 or 3")