Il giorno lunedì 17 febbraio 2020 17:48:07 UTC+1, alberto ha scritto: > Hi, > I would use this script to evaluate fugacity coefficient with PENG-ROBINSON > equation, but I don't understand the correct mode to insert data > > import numpy as np > import matplotlib.pyplot as plt > from scipy.optimize import newton > > R = 8.314e-5 # universal gas constant, m3-bar/K-mol > class Molecule: > """ > Store molecule info here > """ > def __init__(self, name, Tc, Pc, omega): > """ > Pass parameters desribing molecules > """ > #! name > self.name = name > #! Critical temperature (K) > self.Tc = Tc > #! Critical pressure (bar) > self.Pc = Pc > #! Accentric factor > self.omega = omega > > def print_params(self): > """ > Print molecule parameters. > """ > print("""Molecule: %s. > \tCritical Temperature = %.1f K > \tCritical Pressure = %.1f bar. > \tAccentric factor = %f""" % (self.name, self.Tc, self.Pc, > self.omega)) > > def preos(molecule, T, P, plotcubic=True, printresults=True): > """ > Peng-Robinson equation of state (PREOS) > > http://en.wikipedia.org/wiki/Equation_of_state#Peng.E2.80.93Robinson_equation_of_state > :param molecule: Molecule molecule of interest > :param T: float temperature in Kelvin > :param P: float pressure in bar > :param plotcubic: bool plot cubic polynomial in compressibility factor > :param printresults: bool print off properties > Returns a Dict() of molecule properties at this T and P. > """ > # build params in PREOS > Tr = T / molecule.Tc # reduced temperature > a = 0.457235 * R**2 * molecule.Tc**2 / molecule.Pc > b = 0.0777961 * R * molecule.Tc / molecule.Pc > kappa = 0.37464 + 1.54226 * molecule.omega - 0.26992 * molecule.omega**2 > alpha = (1 + kappa * (1 - np.sqrt(Tr)))**2 > > A = a * alpha * P / R**2 / T**2 > B = b * P / R / T > > # build cubic polynomial > def g(z): > """ > Cubic polynomial in z from EOS. This should be zero. > :param z: float compressibility factor > """ > return z**3 - (1 - B) * z**2 + (A - 2*B - 3*B**2) * z - ( > A * B - B**2 - B**3) > > # Solve cubic polynomial for the compressibility factor > z = newton(g, 1.0) # compressibility factor > rho = P / (R * T * z) # density > > # fugacity coefficient comes from an integration > fugacity_coeff = np.exp(z - 1 - np.log(z - B) - A / np.sqrt(8) / B * > np.log( > (z + (1 + np.sqrt(2)) * B) / (z + (1 - np.sqrt(2)) * B))) > > if printresults: > print("""PREOS calculation at > \t T = %.2f K > \t P = %.2f bar""" % (T, P)) > print("\tCompressibility factor : ", z) > print("\tFugacity coefficient: ", fugacity_coeff) > print("\tFugacity at pressure %.3f bar = %.3f bar" % ( > P, fugacity_coeff * P)) > print("\tDensity: %f mol/m3" % rho) > print("\tMolar volume: %f L/mol" % (1.0 / rho * 1000)) > print("\tDensity: %f v STP/v" % (rho * 22.4 / 1000)) > print("\tDensity of ideal gas at same conditions: %f v STP/v" % ( > rho * 22.4/ 1000 * z)) > > if plotcubic: > # Plot the cubic equation to visualize the roots > zz = np.linspace(0, 1.5) # array for plotting > > plt.figure() > plt.plot(zz, g(zz), color='k') > plt.xlabel('Compressibility, $z$') > plt.ylabel('Cubic $g(z)$') > plt.axvline(x=z) > plt.axhline(y=0) > plt.title('Root found @ z = %.2f' % z) > plt.show() > return {"density(mol/m3)": rho, "fugacity_coefficient": fugacity_coeff, > "compressibility_factor": z, "fugacity(bar)": fugacity_coeff * P, > "molar_volume(L/mol)": 1.0 / rho * 1000.0} > > def preos_reverse(molecule, T, f, plotcubic=False, printresults=True): > """ > Reverse Peng-Robinson equation of state (PREOS) to obtain pressure for a > particular fugacity > :param molecule: Molecule molecule of interest > :param T: float temperature in Kelvin > :param f: float fugacity in bar > :param plotcubic: bool plot cubic polynomial in compressibility factor > :param printresults: bool print off properties > Returns a Dict() of molecule properties at this T and f. > """ > # build function to minimize: difference between desired fugacity and > that obtained from preos > def g(P): > """ > :param P: pressure > """ > return (f - preos(molecule, T, P, plotcubic=False, > printresults=False)["fugacity(bar)"]) > > # Solve preos for the pressure > P = newton(g, f) # pressure > > # Obtain remaining parameters > pars = preos(molecule, T, P, plotcubic=plotcubic, > printresults=printresults) > rho = pars["density(mol/m3)"] > fugacity_coeff = pars["fugacity_coefficient"] > z = pars["compressibility_factor"] > > return {"density(mol/m3)": rho, "fugacity_coefficient": fugacity_coeff, > "compressibility_factor": z, "pressure(bar)": P, > "molar_volume(L/mol)": 1.0 / rho * 1000.0} > > # TODO: Implement mixture in object-oriented way as well > def preos_mixture(molecule_a, molecule_b, delta, T, P_total, x, > plotcubic=True, printresults=True): > """ > Peng-Robinson equation of state (PREOS) for a binary mixture > > http://en.wikipedia.org/wiki/Equation_of_state#Peng.E2.80.93Robinson_equation_of_state > :param molecule_a: Molecule molecule 1 of interest > :param molecule_b: Molecule molecule 2 of interest > :param delta: binary interaction parameter between molecules a and b > :param T: float temperature in Kelvin > :param P_total: float total pressure in bar > :param x: array mole fractions > :param plotcubic: bool plot cubic polynomial in compressibility factor > :param printresults: bool print off properties > """ > # build arrays of properties > Tc = np.array([molecule_a.Tc, molecule_b.Tc]) > Pc = np.array([molecule_a.Pc, molecule_b.Pc]) > omega = np.array([molecule_a.omega, molecule_b.omega]) > > # build params in PREOS > Tr = T / Tc # reduced temperature > a0 = 0.457235 * R**2 * Tc**2 / Pc > b = 0.0777961 * R * Tc / Pc > kappa = 0.37464 + 1.54226 * omega - 0.26992 * omega**2 > a = a0 * (1 + kappa * (1 - np.sqrt(Tr)))**2 > > # apply mixing rules > aij = (1.0 - delta) * np.sqrt(a[0] * a[1]) > a_mix = a[0] * x[0]**2 + a[1] * x[1]**2 + 2.0 * x[0] * x[1] * aij > b_mix = x[0] * b[0] + x[1] * b[1] > A = a_mix * P_total / R**2 / T**2 > B = b_mix * P_total / R / T > > # build cubic polynomial > def g(z): > """ > Cubic polynomial in z from EOS. This should be zero. > :param z: float compressibility factor > """ > return z**3 - (1 - B) * z**2 + (A - 2*B - 3*B**2) * z - ( > A * B - B**2 - B**3) > > # Solve cubic polynomial for the compressibility factor > z = newton(g, 1.0) # compressibility factor > rho = P_total / (R * T * z) # density > > Lnfug_0 = -np.log(z - B) + (z - 1.0) * b[0] / b_mix - A / np.sqrt(8) / B > * (2.0 / a_mix * (x[0] * a[0] + x[1] * aij) - b[0] / b_mix) *\ > np.log((z + (1.0 + np.sqrt(2)) * B) / (z + (1.0 - np.sqrt(2)) * B)) > Lnfug_1 = -np.log(z - B) + (z - 1.0) * b[1] / b_mix - A / np.sqrt(8) / B > * (2.0 / a_mix * (x[1] * a[1] + x[0] * aij) - b[1] / b_mix) *\ > np.log((z + (1.0 + np.sqrt(2)) * B) / (z + (1.0 - np.sqrt(2)) * B)) > > # fugacity coefficient comes from an integration > fugacity_coefs = np.exp(np.array([Lnfug_0, Lnfug_1])) > > if printresults: > print("""PREOS calculation at > \t T = %.2f K > \t P, total = %.2f bar""" % (T, P_total)) > print("\tDensity: %f mol/m3" % rho) > print("\tCompressibility factor : %f" % z) > print("Component 0, %s:" % molecule_a.name) > print("\tFugacity coefficient: %f" % fugacity_coefs[0]) > print("\tFugacity: %f bar" % (fugacity_coefs[0] * x[0] * P_total)) > print("Component 1, %s:" % molecule_b.name) > print("\tFugacity coefficient: %f" % fugacity_coefs[1]) > print("\tFugacity: %f bar" % (fugacity_coefs[1] * x[1] * P_total)) > > if plotcubic: > # Plot the cubic equation to visualize the roots > zz = np.linspace(0, 1.5) # array for plotting > > plt.figure() > plt.plot(zz, g(zz), color='k') > plt.xlabel('Compressibility, $z$') > plt.ylabel('Cubic $g(z)$') > plt.axvline(x=z) > plt.axhline(y=0) > plt.title('Root found @ z = %.2f' % z) > plt.show() > return {"density(mol/m3)": rho, "fugacity_coefficients": fugacity_coefs, > "compressibility_factor": z} > > > the readme file says that > > As an example calculation, we consider methane at 65.0 bar and 298.0 K. > Methane has a critical temperature of -82.59 deg. C and a critical pressure > of 45.99 bar. Its accentric factor is 0.011. We first create a methane > molecule object and print its stored parameters: > > import preos > # pass name, Tc, Pc, omega > methane = preos.Molecule("methane", -82.59 + 273.15, 45.99, 0.011) > methane.print_params() > > Could I fix it > > regards > > Alberto
Hi, my code preos in one file preos.py my commands are alberto@HENDRIX ~/PREOS $ python3.5 Python 3.5.2 (default, Oct 8 2019, 13:06:37) [GCC 5.4.0 20160609] on linux Type "help", "copyright", "credits" or "license" for more information. >>> import preos >>> methane = Molecule("methane", -82.59 + 273.15, 45.99, 0.011) Traceback (most recent call last): File "<stdin>", line 1, in <module> NameError: name 'Molecule' is not defined regards Alberto -- https://mail.python.org/mailman/listinfo/python-list