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()