---
tags: Mathematical modeling and optimization
title: 4.1 Properties in the Thermodynamics
---

## 4.1 Properties in the Thermodynamics

Stepping from the mechanical description to involve energy representation, thermodynamics is the only pathway, which describes the state, e.g., pressure and temperature and its relation to energy.

### 4.1.1 Equation of state and the Dalton's law

An equation of state (EoS) is a mathematical model that describes the state properties of a fluid (gas or liquid). It relates pressure $(P)$, volume $(V)$, and temperature $(T)$ of a substance.

#### Ideal Gas Law
The Ideal Gas Law is expressed as:
$$
\displaystyle P\ V= n \ Ru \ T
$$
where:
- $Ru$ is the universal gas constant (8.314 J/(mol·K)), and
- $n$ is the number of moles of the gas, also as $\qquad \displaystyle\frac{\text{mass of gas}}{\text{gas molar weight} (M_w)}$

<br/>

Thus, the Ideal Gas Law also appears in different forms like
$\qquad \qquad \qquad \displaystyle P = \rho \ R \ T\qquad$ with density $\rho$ and $\displaystyle R = \frac{Ru}{M_W}$, or

$\qquad \qquad \qquad \displaystyle P = \frac{Ru \ T}{V_m}\qquad$ with molar volume $\displaystyle V_m = \frac{V}{n}$.

The primary assumptions for ideal gas law are:
- __Gas molecules have no volume__:
The volume occupied by the gas molecules themselves is negligible compared to the volume of the container. This means that the gas molecules are considered point particles with no actual volume.
- __No intermolecular forces__:
There are no attractive or repulsive forces between gas molecules. The only interactions are elastic collisions with each other and the walls of the container.
- __Random motion__:
Gas molecules are in constant, random motion, and they move in straight lines until they collide with another molecule or the walls of the container.
- __Elastic collisions__:
Collisions between gas molecules and between gas molecules and the container walls are perfectly elastic. This means that there is no loss of kinetic energy during collisions; the total kinetic energy of the gas molecules remains constant.
- __Large number of molecules__:
A large number of gas molecules are present, which justifies the statistical averaging of properties like pressure and temperature.
- __Kinetic energy proportional to temperature__:
The average kinetic energy of gas molecules is directly proportional to the absolute temperature (measured in Kelvin). This is expressed mathematically as:
$$\frac{1}{2}mv^2\propto T$$
where $m$ is the mass of a gas molecule, $v$ is its velocity, and $T$ is the absolute temperature.

These assumptions hold reasonably well for many gases under conditions of low pressure and high temperature, where the behavior of real gases approximates that of an ideal gas. However, at high pressures and low temperatures, the deviations from ideal behavior become significant due to the finite volume of gas molecules and the intermolecular forces, necessitating the use of more complex equations of state such as the Van der Waals equation.

#### Van der Waal Equation of state
The Van der Waals equation of state offers a deeper understanding of real gas behavior and phase transitions. It introduces necessary corrections for intermolecular forces and molecular volumes, providing a more accurate description of gas behavior under various conditions. By capturing the critical point and using the Maxwell construction, it helps describe the liquid-vapor phase transition, enhancing our understanding of phase changes in real gases.

The Van der Waals equation is given by:
$$
(P + \frac{a}{V_m^2})(V_m−b)=RT
$$

where $a$ and $b$ stand for _the attraction between particles_ and _the volume excluded by a mole of particles_, respectively.

- Correction for Intermolecular Forces ( $a$ ):
The term $\displaystyle \frac{a}{V_m^2}$ accounts for the attractive forces between gas molecules. These forces reduce the pressure exerted by the gas on the container walls. By adding $\displaystyle \frac{a}{V_m^2}$ to the pressure term, the correct for this reduction shows:
$$
P_{\text{ideal}}=P+\frac{a}{V_m^2}
$$

- Correction for Molecular Volume ( $b$ ):
The term $b$ accounts for the finite volume occupied by gas molecules. The volume available for the gas molecules to move is less than the actual volume of the container. Thus,
$$
V_{\text{ideal}}=V_m−b
$$

The Van der Waals equation can describe the behavior of gases near the critical point, where the gas and liquid phases become indistinguishable. The critical point is characterized by specific values of temperature ($T_c$), pressure ($P_c$), and molar volume ($V_{mc}$).

At the [critical point](https://en.wikipedia.org/wiki/Critical_point_(thermodynamics)), the isotherm on a $P$ vs. $V_m$ plot has a horizontal inflection point
<img src="https://upload.wikimedia.org/wikipedia/commons/thumb/3/3e/Real_Gas_Isotherms.svg/495px-Real_Gas_Isotherms.svg.png" width=45%> <img src="https://useruploads.socratic.org/IdCUlS3SueUcqF19Rs0A_PVT_3D_diagram.png" width=48%>
Thus,

$$
\begin{align}
(\frac{\partial P}{\partial V_m})_T=0 \\
\\
(\frac{\partial^2P}{\partial V_m^2})_T=0
\end{align}
$$

Using these conditions, we can derive expressions for $T_c$, $P_c$, and $V_{mc}$.

Involving the concept of critical point, the van der Waals equation of state is then capable to evaluate the phase state. In every isothermal state (T = const) with a given pressure, for which 3 solutions yields from solving the van der Waals equation, denotes the coexist of liquid and vaper phase.

An example detecting the CO2 state shows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fsolve, newton

# Constants for CO₂

a = 3.59 # L^2 atm / mol^2
b = 0.0427 # L / mol
R = 0.0821 # L atm / K mol


# Function to solve Van der Waals equation for V_m
def van_der_waals_eqn(V_m, T):
    P   = R * T / (V_m - b) - (a / V_m**2)
    return P

# Function to solve Van der Waals equation for V_m
def solv_van_der_waals_eqn(V_m, P, T):
    error = (P + a / V_m**2) * (V_m - b) / (R * T) -1
    return error

# Parameters for the calculation
P = 20 # Pressure in atm
T = 250 # Temperature in K

# Initial guess for molar volume (V_m) in L/mol
guess = (np.linspace(0, 100, 200))*1e-2
# Solve for molar volume using fsolve
Vm = newton(solv_van_der_waals_eqn, guess, args=(P, T), rtol = 1e-10)


Vm_solution = np.unique(Vm.round(6))

Vm_ = []
number_of_solution = 0
for  solution in Vm_solution:
    p_rep = van_der_waals_eqn(solution, T)
    if solution > 0 and p_rep.round(1) == P:
        number_of_solution += 1
        Vm_.append(solution)

if number_of_solution == 1:
    phaseState = "Gas"
elif number_of_solution == 3:
    phaseState = "Liquid and vaper coexist"
else:
    phaseState = "Solution error!"
    

# Calculate compressibility factor Z
Z = (P * Vm_[-1]) / (R * T)

# Output results
print(f"Calculated Molar Volume (V_m): {Vm_} L/mol")
print(f"Compressibility Factor (Z): {Z:.5f}")
print(f"Phase: {phaseState}")

#### Real Gas models and compressibility factor

[Van der Waals equation](https://en.wikipedia.org/wiki/Van_der_Waals_equation) is not the only method describing the real gas behavior, [Redlich-Kwong equation](https://en.wikipedia.org/wiki/Redlich%E2%80%93Kwong_equation_of_state), [Peng-Robinson equation](https://en.wikipedia.org/wiki/Cubic_equations_of_state#Peng%E2%80%93Robinson_equation_of_state) are the conventional method dedicating in this topic, as well as industrial specific description, e.g. [AGA8](https://www.flowsolv.com/examples/density-compressibility-aga8) and [GERG2008](https://pubs.acs.org/doi/10.1021/je300655b)

__The compressibility factor $(\mathbb{Z})$__ is defined as the ratio of the molar volume of a real gas to the molar volume of an ideal gas at the same temperature and pressure. It provides a measure of the deviation of a real gas from ideal gas behavior, where $\mathbb{Z}=1$ for an ideal gas.

### Dalton's Law.

Dalton's Law of Partial Pressures states that in a mixture of non-reacting gases, the total pressure exerted is equal to the sum of the partial pressures of the individual gases. Thus,
$$
P_{\text{total}}=P_1+P_2+ \dots + P_n
$$
where $P_{\text{total}}$ is the total pressure and $P_i$​ is the partial pressure of the $i$-th gas.

The partial pressure of each gas in a mixture can be found using:
$$
P_i=\xi_i\ P_{\text{total}}
$$
where $\xi_i$ is the mole fraction of the $i$-th gas, defined as $\xi_i=\frac{n_i}{n_\text{total}}$

### 4.1.2 Calorific perfect gas assumption

#### Internal energy and enthalpy

__Internal energy ($U$)__ is the total energy contained within a thermodynamic system due to the microscopic motions and interactions of its particles. It includes kinetic energy (due to molecular motion) and potential energy (due to intermolecular forces).

Internal energy is a state function, meaning its value depends only on the state of the system (defined by properties such as temperature, pressure, and volume) and not on how the system arrived at that state.

</br>

__Enthalpy ($H$)__ is a thermodynamic potential that represents the total heat content of a system. It is the sum of the internal energy and the product of pressure and volume:
$$
H = U + P\ V
$$

Enthalpy is particularly useful in processes occurring at constant pressure, such as chemical reactions and phase changes. It simplifies the analysis of heat transfer since at constant pressure, the heat added or removed is directly related to the change in enthalpy.

#### Specific Heat Capacity

Specific heat capacity $(C)$ is defined as the amount of heat required to raise the temperature of a unit mass (or mole) of a substance by one degree Celsius (or one Kelvin). It is an intrinsic property of the material.

The specific heat capacity can be defined at constant pressure $(C_p)$ or constant volume $(C_v)$. The formulation reads:
$$
\begin{align}
C_v = \left( \frac{\partial U}{\partial T} \right)\vert_{v=\text{const}} \\
\\
C_p = \left( \frac{\partial H}{\partial T} \right)\vert_{p=\text{const}} \\
\end{align}
$$
where $U$ and $H$ are the internal energy and enthalphy per unit mass (or mole) of substance.

#### Calorically Perfect Gas
A calorically perfect gas is an idealized gas formulation with constant specific heats. This means that $(C_p)$ and $(C_v)$ do not vary with temperature. While in reality, specific heats can change with temperature, the assumption of constant specific heats is often acceptable for engineering purposes over moderate temperature ranges.

For a calorically perfect gas:
$$
C_p = \text{const}, \qquad \text{and}\qquad C_v = \text{const}
$$

Due to the relation of internal energy and enthalpy, the behavior between $(C_p)$ and $(C_v)$ for an ideal gas is given by:
$$
C_p - C_v = R .
$$

### 4.1.3 First law of thermodynamics

[The __first law of thermodynamics__](https://en.m.wikipedia.org/wiki/First_law_of_thermodynamics) is a statement of the conservation of energy. It asserts that energy cannot be created or destroyed, only transferred or converted from one form to another, which states:
$$
Q = W + \Delta E
$$

Where:

- $\Delta E$ is the change in energy state of the system.
- $Q$ is the heat added to the system.
- $W$ is the work done by the system.

The energy state $(E)$ obtains the internal energy $(U)$ and the kinetic energy $\displaystyle (\frac{1}{2}m u^2)$

In a closed system undergoing a process with negligible kinetic energy, the change in internal energy ($d\ U$) is net of heat added to the system $(\delta \ Q)$ subtracted by the work done by the system $(\delta\ W)$ within a certain time period.
$$
d\ U=\delta\ Q−\delta \ W
$$

For a quasi-static process where the work is done by the system due to pressure $p$ acting through a change in volume $dV$:
$$
\delta\ W=\ p\ dV,
$$
Also, by involving the enthalpy derivative as
$$
d\ H = d\ U + d(p\ V)
$$
and the heat into the system has several representations:
$$
\delta\ Q = d \ U + p\ dV = d \ H− V\ dp
$$

The internal energy and enthalpy for a calorically perfect gas can be expressed as linear functions of temperature due to the constancy of $C_p$ and $C_v$.

The __change in internal energy__ $( \Delta U )$ for a calorically perfect gas by the mole number $(n)$ of a substance yields $\Delta U = n\ C_v\ \Delta T$. This denots the internal energy in the a system with given temperature $(T)$ as
$$
U = n \ C_v\ T
$$

The change in enthalpy $( \Delta H )$ for a calorically perfect gas has the same form.

#### Thermodynamic Processes:
- __Isobaric Process__ (Constant Pressure): the heat added or removed from the system $( Q_p )$ is related to the change in enthalpy: $Q_p = \Delta H = n C_p \Delta T$
- __Isochoric Process__ (Constant Volume): the heat added or removed from the system $( Q_v )$ is related to the change in internal energy: $Q_v = \Delta U = n C_v \Delta T$
- __Isothermal Process__ (Constant Temperature): involving an ideal gas, the change in internal energy is zero $( \Delta U = 0 )$, since internal energy of an ideal gas depends only on temperature. The work done by or on the gas is balanced by the heat added or removed:
$$
Q = W = nRT \ln \left( \frac{V_f}{V_i} \right)
$$
- __Adiabatic Process__ (No Heat Transfer): the first law of thermodynamics simplifies to: $\Delta U = -W$. For a calorically perfect gas undergoing a reversible adiabatic process (also known as an isentropic process), the following relations hold:
$$
\begin{align}
TV^{\gamma-1} =& \text{const} \\
\\
PV^{\gamma} =& \text{const} \\
\end{align}
$$
where $(\gamma)$ is the adiabatic index or heat capacity ratio, defined as
$$
\gamma =\frac{C_p}{C_v}
$$

#### Gas mixing process

Gas mixing process is a common application in this scope. By inserting a certain amount of gas into a vessel, the temperature can change due to the different energy contained. Here is an example:

_Given a 200 liter gas container of 1 bar Air under ambient temperature (18°C). Imagine of an adiabatic system, 5 mole of CO2 by 20°C is injected (isobaric) into the system. What is the final temperature and pressure level in the container?_

import numpy as np
from scipy.optimize import newton

Name_Gas_Composites     =       ["Air", "CO2"]
moleMass_Gas_Composites =       [0.028964, 0.044]        #kg
Cp_Gas_Composites =             [29.19, 36.94]           #J/(mol K)
Cv_Gas_Composites =             [20.85, 28.46]           #J/(mol K)
# https://en.wikipedia.org/wiki/Table_of_specific_heat_capacities
Ru = 8.314                                               #J/(mol K) universal gas constant

vdW_constant_a =        np.array([1.4, 3.640])*1e-1      #[m^6 Pa/mol]
vdW_constant_b =        np.array([0.039, 0.04267])*1e-3  #[m^3/mol]    
#======================     Functions   ==========================#

def Daltons_law(mole_number, temp, method):
    
    # --- Dalton's law --- #
    pressure = 0 
    for it, mole in enumerate(mole_number):
        if method == "van der Waals":
            pressure += van_der_waals_eqn(mole, temp, vdW_constant_a[it], vdW_constant_b[it])
        else:
            pressure += ideal_gas_law(mole, temp)
    
    return pressure
#----------------------------------------------------------------#

def ideal_gas_law(mole, temp):
    
    pressure = mole* Ru * temp /volume
    return pressure

#-----------------------------------------------------------------#
# Function to solve Van der Waals equation for mole number
def solv_moleNumber_igl(mole, pressure, temp):
    
    error = pressure * volume - (mole * Ru * temp)
    return error
#----------------------------------------------------------------#

def van_der_waals_eqn(mole, temp, a, b):
    V_m = volume/mole
    pressure   = Ru * temp / (V_m - b) - (a / V_m**2)
    return pressure

#-----------------------------------------------------------------#
# Function to solve Van der Waals equation for mole number
def solv_moleNumber_vdW(mole, pressure, temp, a, b):
    
    V_m = volume/mole
    error = (pressure + a / V_m**2) * (V_m - b) / (Ru * temp) -1
    return error

#----------------------------------------------------------------#
def calc_internalEnergy(mole_number, temp):
    
    internalEnergy = 0
    for it, mole in enumerate(mole_number):
        internalEnergy += mole*Cv_Gas_Composites[it]*temp
    return internalEnergy

#----------------------------------------------------------------#
def solv_temperature(temp, mole_number, internalEnergy_previous, enthalpy):
    
    internalEnergy_current = calc_internalEnergy(mole_number, temp)
    error = internalEnergy_current - (internalEnergy_previous + enthalpy) 
    return error
#================================================================#


volume = 200*1e-3                                       #[m^3]
pressure = 1.013e5                                      #[Pa]
temperature_initial = 18 + 273.015                      #[K]
mole_number_injected = np.array([0, 5])                #[mole]
temperature_injected = 20+ 273.015                       #[K]
method = "ideal Gas Law"# "van der Waals" #or "Ideal Gas Law"

print(f"---- applying {method} equation -----")
if method == "van der Waals":
    mole_number_temp = newton(solv_moleNumber_vdW, 1e-5, \
                              args=(pressure, temperature_initial, vdW_constant_a[0], vdW_constant_b[0]))
else:
    mole_number_temp = newton(solv_moleNumber_igl, 0.0, args=(pressure, temperature_initial))

                                      
mole_number_initial = np.array([mole_number_temp, 0.0])

internalEnergy_initial =calc_internalEnergy(mole_number_initial, temperature_initial)

enthlapy_injected = sum(mole_number_injected*np.array(Cp_Gas_Composites)*temperature_injected)

mole_number = mole_number_initial + mole_number_injected

temperature_mixed = newton(solv_temperature, temperature_initial, \
                           args=(mole_number, internalEnergy_initial, enthlapy_injected), tol=1e-08 )

    
pressure_mixed = Daltons_law(mole_number, temperature_mixed, method) 
   
print(f"mixed state \n pressure {pressure_mixed/1e5:.3f} bar,\n temperature: {temperature_mixed-273.015:.3f} °C")

Extending the same idea to a serie of gas injection, the tank test application is one of the classical scenarios. __Tank Test__ is a standardized test device/procedure for the airbag inflator design and manufacturing process.

Airbag inflator is the core element in an airbag system, which releases certain amount of gas within the ms time scale. This gas fill up the bag to accomplish the restrain system.

The working principle of tank test is straight forward. The inflator is set and deployed in the closed tank as figure left. During the test, the pressure in the system is measured.

<img src="https://root.hude.com/joomla/images/measure_and_test/produkte/TankTest/th_400_beschriftet.jpg" width = 45%> <img src="https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcS93UY5o6P5KQRnQ9hy8O8VMCLyLQKdPAUssQ&s" width="53%">

_Imagine an inflator deploys with the gas mass flow rate of Helium as 0.8 kg/s, with 30°C for 25ms into the tank. The tank consists of 28.3 liter with initially 23°C filled by air. How would the measured tank curve look like?_

import numpy as np
from scipy.optimize import newton
import matplotlib.pyplot as plt

Name_Gas_Composites     =       ["Air", "He"]
moleMass_Gas_Composites =       [0.028964, 0.006]        #kg
Cp_Gas_Composites =             [29.19, 20.7862]           #J/(mol K)
Cv_Gas_Composites =             [20.85, 12.4717]           #J/(mol K)
# https://en.wikipedia.org/wiki/Table_of_specific_heat_capacities
Ru = 8.314                                               #J/(mol K) universal gas constant

vdW_constant_a =        np.array([1.4, 0.0346])*1e-1      #[m^6 Pa/mol]
vdW_constant_b =        np.array([0.039, 0.0238])*1e-3  #[m^3/mol]    
#======================     Functions   ==========================#

def Daltons_law(mole_number, temp, method):
    
    # --- Dalton's law --- #
    pressure = 0 
    for it, mole in enumerate(mole_number):
        if method == "van der Waals":
            pressure += van_der_waals_eqn(mole, temp, vdW_constant_a[it], vdW_constant_b[it])
        else:
            pressure += ideal_gas_law(mole, temp)
    
    return pressure
#----------------------------------------------------------------#

def ideal_gas_law(mole, temp):
    
    pressure = mole* Ru * temp /volume
    return pressure

#-----------------------------------------------------------------#
# Function to solve Van der Waals equation for mole number
def solv_moleNumber_igl(mole, pressure, temp):
    
    error = pressure * volume - (mole * Ru * temp)
    return error
#----------------------------------------------------------------#

def van_der_waals_eqn(mole, temp, a, b):
    V_m = volume/mole
    pressure   = Ru * temp / (V_m - b) - (a / V_m**2)
    return pressure

#-----------------------------------------------------------------#
# Function to solve Van der Waals equation for mole number
def solv_moleNumber_vdW(mole, pressure, temp, a, b):
    
    V_m = volume/mole
    error = (pressure + a / V_m**2) * (V_m - b) / (Ru * temp) -1
    return error

#----------------------------------------------------------------#
def calc_internalEnergy(mole_number, temp):
    
    internalEnergy = 0
    for it, mole in enumerate(mole_number):
        internalEnergy += mole*Cv_Gas_Composites[it]*temp
    return internalEnergy

#----------------------------------------------------------------#
def solv_temperature(temp, mole_number, internalEnergy_previous, enthalpy):
    
    internalEnergy_current = calc_internalEnergy(mole_number, temp)
    error = internalEnergy_current - (internalEnergy_previous + enthalpy) 
    return error
#================================================================#


volume = 60.0*1e-3                                       #[m^3]
pressure = 1.013e5                                      #[Pa]
temperature_initial = 18 + 273.015                      #[K]
method = "ideal Gas Law"# "van der Waals" #or "Ideal Gas Law"

massFlowRate = np.array([0, 0.8])                                      #[kg/s]
temperature_injected = 30 + 273.015                       #[K]

molarFlowRate = massFlowRate/moleMass_Gas_Composites

timeStep = 1e-4                                         #[s]
timeEnd = 50*1e-3                                       #[s]
numberOfIteration = int(timeEnd/timeStep)

print(f"---- applying {method} equantion -----")
if method == "van der Waals":
    mole_number_temp = newton(solv_moleNumber_vdW, 1e-5, \
                              args=(pressure, temperature_initial, vdW_constant_a[0], vdW_constant_b[0]))
else:
    mole_number_temp = newton(solv_moleNumber_igl, 0.0, args=(pressure, temperature_initial))
                                     
mole_number_initial = np.array([mole_number_temp, 0.0])

timeSequence = []
pressureSequence = []
temperatureSequence = []
for i in range(0, numberOfIteration):  
    currentTime = i* timeStep
    
    internalEnergy_initial =calc_internalEnergy(mole_number_initial, temperature_initial)
    
    if currentTime <= 25*1e-3:
        mole_number_injected = molarFlowRate * timeStep               #[mole]
    else:
        mole_number_injected = np.array([0, 0])

    enthlapy_injected = sum(mole_number_injected*np.array(Cp_Gas_Composites)*temperature_injected)

    mole_number = mole_number_initial + mole_number_injected

    temperature_mixed = newton(solv_temperature, temperature_initial, \
                           args=(mole_number, internalEnergy_initial, enthlapy_injected), tol=1e-08 )
    
    pressure_mixed = Daltons_law(mole_number, temperature_mixed, method) 

    
    timeSequence.append(currentTime*1e3)
    pressureSequence.append((pressure_mixed-1.013e5)/1e3)
    temperatureSequence.append(temperature_mixed)
    
    # update for the next iteration
    mole_number_initial = mole_number
    temperature_initial = temperature_mixed

# Plot results
fig, (ax1, ax2) = plt.subplots(2)
fig.suptitle('Tank test results')
ax1.plot(timeSequence, pressureSequence)
ax1.grid(True)
ax1.set_ylabel("pressure [kPa]")

ax2.plot(timeSequence, temperatureSequence)
ax2.grid(True)
ax2.set_xlabel("time [ms]")
ax2.set_ylabel("temperature [°C]")