K calcu.py

import os
import math as m
import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from decimal import Decimal
from scipy.optimize import newton
import itertools


import preos
from preos import Molecule


fugcoeff_peng = []
fugcoeff_van = []
sub_dict = {
    "O2" : [154.6, 50.46, 0.021],
    "SO2" : [430.8, 78.83, 0.251],
    "SO3" : [491.0, 82.07, 0.41],
    "fugcoeff_peng" : fugcoeff_peng,
    "fugcoeff_van" : fugcoeff_van,
}
result_dict = {}
#for P,Sub in itertools.product(range(1,21,1),['O2','SO2','SO3']):
 #   sub_coeff = sub_dict[Sub]
  #  O2 = preos.Molecule(Sub, sub_coeff[0], sub_coeff[1], sub_coeff[2])
   # props = preos.preos(O2, 690 + 273, P*1.01325 , plotcubic=False, printresults=False)
    #fugcoeff_peng.append(props["fugacity_coefficient_peng"])
    #fugcoeff_van.append(props["fugacity_coefficient_van"])
    #for type in ["peng",'van']:
     #   result_dict[f"{Sub}_{type}"] = sub_dict[f"fugcoeff_{type}"]



def solve_for_Kp(P,Kf=8.373,Ky=1):
    Fug_Coeff_List = []
    component_P = find_P_of_each_sub(P,Kf,Ky)
    print(component_P)
    for P_comp, Comp in zip(component_P, ["SO2","O2","SO3"]):
        #print(P_comp, Comp)
        sub_coeff = sub_dict[Comp]
        #print(sub_coeff)
        O2 = preos.Molecule(Comp, sub_coeff[0], sub_coeff[1], sub_coeff[2])
        props = preos.preos(O2, 690 + 273, P_comp, plotcubic=False, printresults=False)
        Fug_Coeff_List.append(props["fugacity_coefficient_peng"])
    
    Ky_new = (Fug_Coeff_List[2]**2) / ((Fug_Coeff_List[0]**2)*(Fug_Coeff_List[1]))
    if (Ky_new/Ky - 1) == 0:
        return Kf/Ky
    else:
        return solve_for_Kp(P,Kf=8.373, Ky=Ky_new)




def find_P_of_each_sub(P,Kf,Ky):
    dN = -1
    Kp = Kf/Ky # Kf to Kp
    y = sp.symbols('y')
    equation = sp.Eq((((y**2)*((1.5/(3-y/2))**dN))/(((2-y)**2)*(1-y/2))), Kp)
    ans = sp.solve(equation, y)
    for i in range(len(ans)):
        if ans[i].is_real and ans[i] < 2:
            SO3 = ans[i]
    SO2 = 2 - SO3
    N2 = 3.762
    O2 = 1 - SO3/2
    Xso3 = SO3/(O2+SO2+SO3+N2)
    Xso2 = SO2/(O2+SO2+SO3+N2)
    Xo2 = O2/(O2+SO2+SO3+N2)
    Xn2 = N2/(O2+SO2+SO3+N2)
    Pso2 = round(Xso2*P*1.01325,4)
    Po2 = round(Xo2*P*1.01325,4)
    Pso3 = round(Xso3*P*1.01325,4)
    return [float(Pso2), float(Po2), float(Pso3)]
    

def KCal_Call(P):
    Kp_calculated = []
    for P in Prange:
        Kp_called = solve_for_Kp(P)
        Kp_calculated.append(Kp_called)
    print(Kp_calculated)
    return Kp_calculated


#plt.plot(a)
##plt.show()

#
P = 1
Ky = 1
dN = -1

Kf = 8.373#699417047245
Kp = Kf//Ky
Kp_Value = Kp
y = sp.symbols('y')
equation = sp.Eq((((y**2)*((1.5/(3-y/2))**dN))/(((2-y)**2)*(1-y/2))), Kp_Value)
ans = sp.solve(equation, y)
for i in range(len(ans)):
    if ans[i].is_real and ans[i] < 2:
        SO3 = ans[i]
SO2 = 2 - SO3
N2 = 3.762
O2 = 1 - SO3/2
Xso3 = SO3/(O2+SO2+SO3+N2)
Xso2 = SO2/(O2+SO2+SO3+N2)
Xo2 = O2/(O2+SO2+SO3+N2)
Xn2 = N2/(O2+SO2+SO3+N2)
Pso2 = Xo2*P
Po2 = Xso2*P
Pso3 = Xso3*P
#print(SO2, O2, SO3)
#print(Xso2,Xo2,Xso3,Xn2)
#print(Pso2,Po2,Pso3)

#a = Gamma(range(1,21,1), [690],"atm","C")

#print("Fugacity: ", a)


#







T2 = 690 + 273
R = 8.314
dG2 = T2*(-((3.3256*(10**-8)*(T2**3) + 205.485*T2 - 747.844)/(T2**2)) - 0.00900406*m.log(T2) + 0.256796)*1000
print("\nDelta G at T2 =",dG2)
Kp_T2 = m.exp(-dG2/(R*T2))
print("Kp = ",Kp_T2,"\n")











Kp = Kp_T2
dN = -1
R = 8.314
T1 = 298


T_List = np.arange(300,751,25)
T_List_C = np.arange(300,751,50)
dG1 = -141.74*10**3
dH1 = -197.78*10**3
Kp_L = []
dG2_list = []
for T2 in T_List:
    T2 += 273
    dG2 = T2*(-((3.3256*(10**-8)*(T2**3) + 205.485*T2 - 747.844)/(T2**2)) - 0.00900406*m.log(T2) + 0.256796)*1000
    Kp_T2 = m.exp(-dG2/(R*T2))
    dG2_list.append(dG2)
    Kp_L.append(Kp_T2)
conversions = []
solution=[]
y_list = []
for Kp_Value in Kp_L:
    y = sp.symbols('y')
    equation = sp.Eq((((y**2)*((1.5/(3-y/2))**dN))/(((2-y)**2)*(1-y/2))), Kp_Value)
    solution.append(sp.solve(equation, y))
y_list = solution
for Solutions in solution:
    for i in range(len(Solutions)):
        if Solutions[i].is_real and Solutions[i] < 2:
            conversions.append(float(Solutions[i]/2))
#print(conversions,"-Conv\n")
#print(dG2_list,"-dG2\n")
#print(Kp_L,"-Kp\n")
#print(max(conversions))





    
dN = -1
List_of_Kp = Kp_calculated
List_of_Press = List_of_Kp
#print(Kp_L)
#List_of_Kp = [8.373, 8.373, 8.373, 8.373, 8.373, 8.373, 8.373, 8.373, 8.544, 8.632, 8.7218]
#List_of_Press = [1, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20]
Kx = [Kp/(m.pow(P,dN)) for Kp,P in zip(List_of_Kp,List_of_Press)]
Kc = [Kp/m.pow(R*T2,dN) for Kp in List_of_Kp]
Kp = List_of_Kp

plt.xlabel("P (atm)",fontsize=12)
plt.ylabel("Kp",fontsize=12)
plt.title("Kp",fontsize=12)
plt.xticks(List_of_Press,fontsize=12)
plt.yticks(Kp,fontsize=12)
plt.plot(List_of_Press, np.array(Kp))
plt.grid()
plt.show()

plt.xlabel("P (atm)",fontsize=12)
plt.ylabel("Kx",fontsize=12)
plt.title("Kx",fontsize=12)
plt.xticks(List_of_Press,fontsize=12)
plt.yticks(Kx,fontsize=12)
plt.plot(List_of_Press, np.array(Kx))
plt.grid()
#plt.show()

plt.xlabel("P (atm)",fontsize=12)
plt.ylabel("Kc",fontsize=12)
plt.title("Kc",fontsize=12)
plt.xticks(List_of_Press,fontsize=12)
plt.yticks(Kc,fontsize=12)
plt.plot(List_of_Press, np.array(Kc))
plt.grid()
#plt.show()



#print(conversions)
plt.xlabel("T (K)",fontsize=12)
plt.ylabel("Conversion of SO2",fontsize=12)
plt.title("Conversion of SO2",fontsize=12)
plt.xticks(T_List,fontsize=12)
plt.yticks(np.linspace(max(conversions),min(conversions),10), fontsize=12)
#plt.yticks(conversions, fontsize=12)
plt.plot(T_List, np.array(conversions))
plt.grid()
#plt.show()


#figure, axis = plt.subplots(1,2,squeeze=False)

#axis[0,0].plot(np.arange(1,21,1),Kx)
#axis[0,0].set_title("Kx")
#axis[0,0].set_xlabel("P (atm)")
#axis[0,0].set_ylabel("Kx")
#axis[0,0].set_xticks(range(1,21,1))
#axis[0,0].set_yticks(Kx)
#axis[0,0].grid()


#axis[0,1].plot(T_List, conversions)
#axis[0,1].set_title("Conversions of SO2")
#axis[0,1].set_xlabel("T (K)")
#axis[0,1].set_ylabel("Conversions of SO2")
#axis[0,1].set_xticks(T_List_C)
#axis[0,1].grid()

##plt.show()