CalebBell / thermo

Thermodynamics and Phase Equilibrium component of Chemical Engineering Design Library (ChEDL)
MIT License
594 stars 114 forks source link

Phase diagram for pure components (3 phases) #103

Closed OdhranOOC closed 2 years ago

OdhranOOC commented 2 years ago

Hi Caleb,

I was hoping to generate TP phase diagrams (solid, liquid, gas) for some pure components (CO2,N2,O2). I can do this for two phases (liquid and gas) using Tbubble/Tdew (simply calculating the boiling point for a range of pressures and plotting), however I am unsure of how to extend this to the solid phase. Is this something thermo can do?

I have checked through the previous issues/questions on this page and have not found any replicas, however I did find one other post that deals with mixture phases which might be useful. Here is the code from this question:

from thermo import ChemicalConstantsPackage, CEOSGas, CEOSLiquid, PRMIX, FlashVL
from thermo.interaction_parameters import IPDB
constants, properties = ChemicalConstantsPackage.from_IDs(['co2', 'nitrogen', 'oxygen'])
kijs = IPDB.get_ip_asymmetric_matrix('ChemSep PR', constants.CASs, 'kij')
kijs

eos_kwargs = {'Pcs': constants.Pcs, 'Tcs': constants.Tcs, 'omegas': constants.omegas, 'kijs': kijs}
gas = CEOSGas(PRMIX, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
liquid = CEOSLiquid(PRMIX, eos_kwargs=eos_kwargs, HeatCapacityGases=properties.HeatCapacityGases)
flasher = FlashVL(constants, properties, liquid=liquid, gas=gas)
zs = [0.44, 0.4, 0.16]
flasher.debug_PT(zs=zs, Tmin=100, Tmax=300, Pmin=1e5, Pmax=1e6, pts=50)

Is it possible to adapt this code so that 3 phase PT diagrams can be made for some pure components?

Also for reference here is the code I used to generate the TP curves for liquid/gas phases:

from thermo.chemical import Mixture
import numpy as np
import matplotlib.pyplot as plt
import pprint as pp

bar2pa = 101367
pressures = np.linspace(50000,3000000, num=1000, endpoint=True)
CO2Tbub = np.empty(pressures.__len__())
N2Tbub = np.empty(pressures.__len__())
O2Tbub = np.empty(pressures.__len__())
mixTbub = np.empty(pressures.__len__())
mixTdew = np.empty(pressures.__len__())

for i, pres in enumerate(pressures):
    try:
        co2 = Mixture(['CO2'], P = pres)
        CO2Tbub[i] = co2.Tdew
        N2 = Mixture(['nitrogen'], P=pres)
        N2Tbub[i] = N2.Tdew
        O2 = Mixture(['oxygen'], P=pres)
        O2Tbub[i] = O2.Tdew
        air_co2 = Mixture(['CO2', 'nitrogen', 'oxygen'], ws = [0.44696, 0.4, 1-0.4-0.44696], P = pres)
        mixTbub[i] = air_co2.Tbubble
        mixTdew[i] = air_co2.Tdew
    except:
        pass

plt.figure()
plt.plot(CO2Tbub, pressures, label = 'CO2')
plt.plot(N2Tbub, pressures, label = 'N2')
plt.plot(O2Tbub, pressures, label = 'O2')
plt.plot(mixTbub, pressures, label = 'T bubble Mixture')
plt.plot(mixTdew, pressures, label = 'T bubble Mixture')
plt.legend()
plt.xlabel('Temperature [k]')
plt.ylabel('Psat (MPa)')
plt.legend()
plt.show()

Edit: Also as an unrelated side note, I recently found your Jupiter notebooks with example uses of thermo. Really useful stuff. Thank you for these

CalebBell commented 2 years ago

Hi Odhran, I started to work on solids using a simple sublimation pressure model at one point but it is very much not complete. Thermo doesn't have this functionality today. You can always take VL results from thermo and add the solid yourself in a tool like Excel. Caleb

OdhranOOC commented 2 years ago

You can always take VL results from thermo and add the solid yourself in a tool like Excel.

Thanks Caleb. Are you saying to use excels curve fitting functionality combined with VL data to produce a polynomial sublimation curve?

CalebBell commented 2 years ago

Hi Odhran, Yes, that's right. There is no sublimation data in Thermo at this time. If you find it yourself in the literature, you can combine thermo's VL results manually with the sublimation curve. Sincerely, Caleb

OdhranOOC commented 2 years ago

Perfect thank you 👍