Untitled

unknown
plain_text
2 years ago
3.6 kB
3
Indexable
Never
# -*- coding: utf-8 -*-
"""
Created on Thu Nov 11 15:23:11 2021

@author: carlo
"""

import thermo
import matplotlib.pyplot as plt
import numpy as np

#Set the properties for our substance
methane = np.array([190.6, 4.6e6, .008])

#Code for problem 1
def problem_1():
#initial conditions
P_0, T_0, R = .48e6, 140, 8.31446
param = thermo.PR_param(P_0, T_0, methane)
A, B = param[3], param[4]
z_values=thermo.PR_Z(A,B)
#obtain z-values, intialize our correct z
zl, zh, z  = z_values[0] , z_values[2], 0
#determine the smallest z
if thermo.PR_log_phi(A,B, zl) < thermo.PR_log_phi(A,B, zh):
z = zl
else:
z = zh
print("Z of Methane = " + str(z))
#Calculate the molar volume and enthalpy function
mVol = (z * R * T_0) / P_0
h_d = thermo.PR_H(P_0, T_0, param, z)
print("Molar volume = " + str(mVol))
print("Enthalpy departure = " + str(h_d))

#Listed below are my functions for problem 2
#calculates difference between the fugacities
def df(guess, t):
#find z values
params = thermo.PR_param(guess, t, methane)
A, B = params[3], params[4]
z_values = thermo.PR_Z(A,B)
zl, zh  = z_values[0] , z_values[2]
#calculate the fugacities
fl = np.exp(thermo.PR_log_phi(A,B, zl)) * guess
fh = np.exp(thermo.PR_log_phi(A,B, zh)) * guess
#take the absolute value of the ratio of the fugacities minus one
df = abs((fl/fh)-1)
fugacity_ratio = fl/fh
#returns this quanitity and the ratio of the fugacities
return df, fugacity_ratio

#returns the approximate vapor pressure
def vp(guess, t):
iterations = 1
#checks if our error is within a small value
while abs(df(guess,t)[0]) > 0.0001:
guess *= df(guess,t)[1]
iterations += 1
#returns the vapor pressure in MPa and the number of iterations
return guess / 1e6 , iterations

#a simple function to print vapor pressure data
def print_vp(guess, t):
print("Temperature:", t)
print("Approximated Vapor Pressure:",  vp(guess,t)[0])
print("Iterations:", vp(guess,t)[1])
print("")

#prints out the vapor pressures for the 5 given temperatures
def problem_2():
print_vp(0.05e6,100)
print_vp(.1e6,120)
print_vp(.3e6,140)
print_vp(.5e6,160)
print_vp(2.5e6,180)

#makes two semilog plots: one w a t axis and another with a t/1000 axis
def problem_3():
#initialize the x and y arrays for each plot
temps = np.linspace(100,180, 5)
real_pressures = np.array([0.034931,0.199173,0.686465,1.65732,3.28685])
approx_pressures = np.array([vp(0.05e6,100)[0],
vp(.1e6,120)[0],vp(.3e6,140)[0],vp(.5e6,160)[0],vp(2.5e6,180)[0]])

#creates the semilog plot of vapor pressure with respect to t
plt.figure(figsize=(12,4))
plt.subplot(1,2,1)
plt.semilogy(temps,approx_pressures, 'r' , label = "Approximation")
plt.semilogy(temps,real_pressures, label = "Real")
plt.title("Real vs Approximated Vapor Pressures")
plt.xlabel("Temperature (K)")
plt.ylabel("Vapor Pressure (MPa)")
plt.legend()

#creates the semilog plot of vapor pressure with respect to t/1000
plt.subplot(1,2,2)
plt.semilogy(1000 / temps,approx_pressures, 'r' , label = "Approximation")
plt.semilogy(1000/ temps,real_pressures, label = "Real")
plt.title("Real vs Approximated Vapor Pressures")
plt.xlabel("Temperature (K)")
plt.ylabel("Vapor Pressure (MPa)")
plt.legend()

problem_1()
problem_2()
problem_3()