Solving thermodynamics problems
Contents
Solving thermodynamics problems#
This module introduces how to solve thermodynamics problems in Python using Cantera and Pint. It will also briefly show how CoolProp could also be used, if you need access to a wider range of fluids than Cantera currently supports.
For help running these examples or setting up your own problems, see the module on “Setting up your computing environment”.
First, we need to import the necessary libraries:
# Numpy adds some useful numerical types and functions
import numpy as np
# Cantera will handle thermodynamic properties
import cantera as ct
# Pint gives us some helpful unit conversion
from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity # We will use this to construct quantities (value + unit)
Quick Pint intro#
Pint is a useful tool for parsing unit expressions, converting between different units, and working with expressions involving quantities (values + units).
For example, we can read in expressions that give temperature and pressure in US customary units (e.g., imperial units) and then convert those to SI units:
Pint interprets nearly all properties either written out (e.g., kelvin, meters) or as abbreviations (e.g., K, m)—except for temperatures. In that case, we need to use kelvin or K, celsius or degC, fahrenheits or degF, and rankine or degR.
temperature = Q_(100, 'degF')
pressure = Q_('10 psi')
new_temperature = temperature.to('K')
new_pressure = pressure.to('Pa')
print(f'Temperature in SI: {new_temperature: .2f}')
print(f'Pressure in SI: {new_pressure: .2f}')
Temperature in SI: 310.93 kelvin
Pressure in SI: 68947.57 pascal
Another way to specify units is to multiply the value by an object with the unit (as a member of the ureg object):
distance = 10 * ureg.meter
print(distance)
10 meter
Pint also handles mathematical operations between physical quantities:
distance1 = Q_(1, 'm')
distance2 = Q_(10, 'cm')
print(distance1 + distance2)
1.1 meter
mass = Q_(1, 'kg')
velocity = Q_(2, 'm/s')
kinetic_energy = 0.5 * mass * velocity**2
print(kinetic_energy)
print(kinetic_energy.to('joule'))
2.0 kilogram * meter ** 2 / second ** 2
2.0 joule
Using Cantera for thermodynamic properties#
Cantera is a suite of tools for solving problems involving chemical kinetics, thermodynamics, and transport processes. We’ll use it here primarily for evaluating thermodynamic properties of fluids.
Cantera comes with built-in liquid/vapor equations of state for multiple fluids:
water:
ct.Water()nitrogen:
ct.Nitrogen()methane:
ct.Methane()hydrogen (H2):
ct.Hydrogen()oxygen (O2):
ct.Oxygen()carbon dioxide (CO2):
ct.CarbonDioxide()n-heptane (C7H16):
ct.Heptane()R134a:
ct.Hfc134a()
Let’s take a look at water, and evaluate its properties at different conditions:
# Create an object to hold the thermodynamic state of water
f = ct.Water()
# Fix the thermodynamic state by specifying temperature and
# specific volume, in SI units (K and m^3/kg)
f.TV = 673.15, 1e-2
# Evaluating the object provides a summary of its properties at this state
f()
water:
temperature 673.15 K
pressure 1.9936e+07 Pa
density 100 kg/m^3
mean mol. weight 18.016 kg/kmol
vapor fraction 1
phase of matter supercritical
1 kg 1 kmol
--------------- ---------------
enthalpy -1.315e+07 -2.3692e+08 J
internal energy -1.335e+07 -2.4051e+08 J
entropy 9078.4 1.6356e+05 J/K
Gibbs function -1.9261e+07 -3.4701e+08 J
heat capacity c_p 6284.6 1.1322e+05 J/K
heat capacity c_v 2694 48535 J/K
We can do the same thing to get the properties of R134a; this time, let’s use Pint to interpret the temperature and pressure as given in units of °C and bar, then convert to SI as needed by Cantera:
f = ct.Hfc134a()
# Specify temperature and pressure
temp = Q_(50, 'degC')
pres = Q_(3, 'bar')
# Fix the thermodynamic state using temperature and pressure
# Remember, Cantera requires SI units
f.TP = temp.to('K').magnitude, pres.to('Pa').magnitude
# see overview of properties
f()
HFC-134a:
temperature 323.15 K
pressure 3e+05 Pa
density 11.944 kg/m^3
mean mol. weight 102.03 kg/kmol
vapor fraction 1
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy 2.4319e+05 2.4813e+07 J
internal energy 2.1807e+05 2.225e+07 J
entropy 1607.4 1.64e+05 J/K
Gibbs function -2.7623e+05 -2.8184e+07 J
heat capacity c_p 917.05 93568 J/K
heat capacity c_v 813.22 82975 J/K
Once we have the thermodynamic state fixed for a fluid, we can easily obtain any properties we need:
# get internal energy in SI units (mass basis)
print(f'Internal energy: {f.u: .2f} J/kg')
# enthalpy
print(f'Enthalpy: {f.h: .2f} J/kg')
# constant pressure specific heat
print(f'c_p: {f.cp: .2f} J')
# constant volume specific heat
print(f'c_v: {f.cv: .2f} J')
# density
print(f'density: {f.density: .2f} kg/m^3')
# critical properties
print(f'critical temperature: {f.critical_temperature: .2f} K')
print(f'critical pressure: {f.critical_pressure: .2f} Pa')
print(f'critical density: {f.critical_density: .2f} kg/m^3')
Internal energy: 218070.47 J/kg
Enthalpy: 243187.37 J/kg
c_p: 917.05 J
c_v: 813.22 J
density: 11.94 kg/m^3
critical temperature: 374.21 K
critical pressure: 4059280.00 Pa
critical density: 511.95 kg/m^3
We can fix the thermodynamic state by any combination of two properties:
# specify pressure and internal enegy
pres = Q_(250e3, 'Pa')
internal_energy = Q_(300e3, 'J/kg')
f.UP = internal_energy.to('J/kg').magnitude, pres.to('Pa').magnitude
# get specific volume
print(f'Specific volume: {f.v: 0.4f} m^3/kg')
Specific volume: 0.1332 m^3/kg
Convenience function to get values in SI units#
It can be a pain (and requires more writing) to constantly convert to SI units and extract the magnitude, when using Pint quantities in combination with Cantera and CoolProp.
We can define a handy convenience function to make this easier:
# for convenience:
def to_si(quant):
'''Converts a Pint Quantity to magnitude at base SI units.
'''
return quant.to_base_units().magnitude
You’ll need to define this function in any Python code or Jupyter notebook when you want to use it, but then it makes things much easier. For example:
temperature = Q_(30, 'degC')
pressure = Q_(1, 'atm')
f = ct.Water()
f.TP = to_si(temperature), to_si(pressure)
f()
water:
temperature 303.15 K
pressure 1.0132e+05 Pa
density 995.73 kg/m^3
mean mol. weight 18.016 kg/kmol
vapor fraction 0
phase of matter liquid
1 kg 1 kmol
--------------- ---------------
enthalpy -1.5845e+07 -2.8546e+08 J
internal energy -1.5845e+07 -2.8546e+08 J
entropy 3956.8 71286 J/K
Gibbs function -1.7044e+07 -3.0707e+08 J
heat capacity c_p 4178.6 75282 J/K
heat capacity c_v 4116.1 74156 J/K
Compare ideal and real gas behavior for nitrogen#
Using Cantera’s built-in pure fluid model for nitrogen, and its ideal gas model for air (which includes nitrogen), we can see an example of the error involved with treating fluids as ideal gases:
temp1 = Q_(500, 'K')
temp2 = Q_(1270, 'K')
pres = Q_(500, 'kPa')
# Real gas
air_real1 = ct.Nitrogen()
air_real1.TP = temp1.to('K').magnitude, pres.to('Pa').magnitude
air_real2 = ct.Nitrogen()
air_real2.TP = temp2.to('K').magnitude, pres.to('Pa').magnitude
delta_u_real = Q_(air_real2.u - air_real1.u, 'J/kg')
print(f'Δu for real gas: {delta_u_real: .2f}')
# ideal gas
air_ideal1 = ct.Solution('air.cti')
air_ideal1.TPX = temp1.to('K').magnitude, pres.to('Pa').magnitude, 'N2:1.0'
air_ideal2 = ct.Solution('air.cti')
air_ideal2.TPX = temp2.to('K').magnitude, pres.to('Pa').magnitude, 'N2:1.0'
delta_u_ideal = Q_(air_ideal2.u - air_ideal1.u, 'J/kg')
print(f'Δu for ideal gas: {delta_u_ideal: .2f}')
diff = 100*np.abs(delta_u_ideal-delta_u_real)/delta_u_real
print(f'difference: {diff.magnitude: .2f}%')
Δu for real gas: 647294.70 joule / kilogram
Δu for ideal gas: 648460.39 joule / kilogram
difference: 0.18%
/tmp/ipykernel_4736/2667467038.py:16: DeprecationWarning: XML_Node::build:
The CTI and XML input file formats are deprecated and will be removed in
Cantera 3.0. Use 'cti2yaml.py' or 'ctml2yaml.py' to convert CTI or XML input
files to the YAML format. See https://cantera.org/tutorials/legacy2yaml.html
for more information.
air_ideal1 = ct.Solution('air.cti')
Fortunately, at these conditions, nitrogen behaves much like an ideal gas, and so our calculations of the internal energy difference between two states are quite close.
Modeling air#
Air is a pseudo-pure fluid; in reality it is a multicomponent gas mixture of mostly oxygen and nitrogen, with trace amounts of arfon, carbon dioxide, and water (depending on the humidity). However, since the composition does not change under most processes, we can sometimes treat it like a pure fluid.
Using Cantera#
To solve problems with air, we can use an ideal gas model in Cantera, or the pseudo-pure fluid model in CoolProp (see the CoolProp module for more details).
In Cantera, we create a Solution object using the built-in air.cti model file, then specify the state, including composition. We can use mole fractions (X) or mass fractions (Y). The commonly accepted composition of air is 1 mole of oxygen to 3.76 moles of nitrogen, and we can specify this molar composition using the string 'O2:1.0, N2:3.76':
air = ct.Solution('air.cti')
air.TPX = 300.0, 101325.0, 'O2:1.0, N2:3.76'
air()
print(f'Density: {air.density: .3f} kg/m^3')
air:
temperature 300 K
pressure 1.0133e+05 Pa
density 1.172 kg/m^3
mean mol. weight 28.851 kg/kmol
phase of matter gas
1 kg 1 kmol
--------------- ---------------
enthalpy 1907.6 55035 J
internal energy -84548 -2.4393e+06 J
entropy 6891.7 1.9883e+05 J/K
Gibbs function -2.0656e+06 -5.9594e+07 J
heat capacity c_p 1010.1 29141 J/K
heat capacity c_v 721.87 20827 J/K
mass frac. Y mole frac. X chem. pot. / RT
--------------- --------------- ---------------
O2 0.233 0.21008 -26.234
N2 0.767 0.78992 -23.269
[ +6 minor] 0 0
Density: 1.172 kg/m^3
A Solution object can be used to specify and compute thermodynamic, chemical kinetic, and/or transport properties of a mixture; the thermodynamics come from the ThermoPhase class.
We can set the state by using a combination of two intensive properties, from the list [DP, HP, SP, SV, TD, TP, UV], along with the mole or mass fractions. Here, D is density, H is enthalpy, P is pressure, S is entropy, V is specific volume, and T is temperature; all should be in SI units, either in a mass or molar basis.
For example, we could specify the temperature, pressure, and air mole fractions using
air.TPX = 300, 101325, 'O2:1.0, N2:3.76'
or we could specify via density and pressure:
air.DPX = 1.225, 101325, 'O2:1.0, N2:3.76'
(Remember that density and specific volume are inversely related, so if we have one we always have the other.)
Using CoolProp#
We can also use CoolProp to model air, using its psuedo-pure fluid model. In this case, we do not need to specify the composition.
from CoolProp.CoolProp import PropsSI
density = PropsSI('D', 'T', 300.0, 'P', 101325, 'air')
print(f'Density: {density: .3f} kg/m^3')
Density: 1.177 kg/m^3
Hmm, the density given by CoolProp differs slightly from that given by Cantera. The difference is less than 1%, but where could that be coming from?
Well, our specified moles for the composition are only approximate, and also missing a small amount of argon. Let’s use a slightly more accurate mixture (though this is still missing traces of carbon dioxide, neon, helium, etc., that are present at less than 0.01% by mole):
air = ct.Solution('air.cti')
air.TPX = 300.0, 101325.0, 'N2:0.78084, O2:0.20947, Ar:0.00934'
print(f'Density: {air.density: .3f} kg/m^3')
Density: 1.176 kg/m^3