Solving the First Law of Thermodynamics using Python

The first law of thermodynamics is related to energy conservation if energy in one form disappears then the energy will appear in some other form. In thermodynamics we are primarily concerned about heat, work and internal energy. Heat and work are forms of energy called "energy in transit", making them path functions that cannot be stored by a system, whereas internal energy is a property of the system that can be stored.

For a closed system, the first law is written as

$$\mathrm{\Sigma Q=\Sigma W}$$

This is only valid for a closed system undergoing a cycle. For a process it can be written as

$$\mathrm{Q_{1-2}=\triangle U+W_{1-2}}$$

Before we can start modeling these equations, we need to understand how to model work and heat for different thermodynamic processes.

Thermodynamic Work and Processes

If the only impact on the system's external environment can be reduced to weight increase, then the system is said to have performed thermodynamic work. Work is a path function since it depends on the path taken between two states. We focus on non-dissipative displacement work, represented by the area under the p-v graph for a closed system.

For a general polytropic process represented by $pv^{n}$=Constant, the work done depends on the index n. The table below shows different processes based on index n

Work interaction is derived by integrating the process from initial to final state

$$\mathrm{W=\int ^{2}_{1}pdv}$$

The work interactions for different processes are shown below

Heat Transfer

Heat transfer is energy transition between a system and surroundings due to temperature difference. Heat entering the system is positive, while heat leaving is negative. When no heat exchange occurs, it's an adiabatic process.

For sensible heat transfer

$$\mathrm{Q=mc\Delta T}$$

For latent heat during phase change

$$Q=mL$$

Python Functions for Thermodynamic Calculations

The following Python functions calculate work and heat interactions for different processes

import numpy as np
import matplotlib.pyplot as plt

# Constant Volume Process
def w_cv(p1, p2, v):
    """
    Function for constant volume process
    Input: Pressure limits p1 & p2 (in kPa) and volume (m^3)
    Output: Work output and p-v plot
    """
    plt.figure(figsize=(8, 5))
    p = np.linspace(min(p1, p2), max(p1, p2), 30)
    v_array = np.zeros(30) + v
    
    plt.plot(v_array, p, 'k-', linewidth=3)
    plt.xlabel('Volume (m³)')
    plt.ylabel('Pressure (kPa)')
    plt.title('Constant Volume Process')
    plt.grid(True, alpha=0.3)
    plt.show()
    return 0

# Constant Pressure Process
def w_cp(p, v1, v2):
    """
    Function for constant pressure process
    Input: Volume limits v1 & v2 (m^3) and pressure (kPa)
    Output: Work output and p-v plot
    """
    work = p * (v2 - v1)
    plt.figure(figsize=(8, 5))
    v = np.linspace(min(v1, v2), max(v1, v2), 30)
    p_array = np.zeros(30) + p
    
    plt.plot(v, p_array, 'k-', linewidth=3)
    plt.fill_between(v, p_array, color='red', alpha=0.3)
    plt.xlabel('Volume (m³)')
    plt.ylabel('Pressure (kPa)')
    plt.title('Constant Pressure Process')
    plt.grid(True, alpha=0.3)
    plt.show()
    return round(work, 3)

# Isothermal Process
def w_ct(p1, v1, v2):
    """
    Function for constant temperature process
    Input: Initial pressure (kPa) and volume (m^3) and final volume (m^3)
    Output: Work output and p-v plot
    """
    work = p1 * v1 * np.log(v2 / v1)
    plt.figure(figsize=(8, 5))
    
    v = np.linspace(min(v1, v2), max(v1, v2), 30)
    c = p1 * v1
    p = c / v
    
    plt.plot(v, p, 'k-', linewidth=3)
    plt.fill_between(v, p, color='red', alpha=0.3)
    plt.xlabel('Volume (m³)')
    plt.ylabel('Pressure (kPa)')
    plt.title('Isothermal Process')
    plt.grid(True, alpha=0.3)
    plt.show()
    return round(work, 3)

# Polytropic Process
def w_pol(p1, p2, v1, v2, n):
    """
    Function for polytropic process
    Input: p1 & p2 (kPa), v1 & v2 (m^3), and n
    If adiabatic process, use n=1.4
    Output: Work output and p-v plot
    """
    work = (p1 * v1 - p2 * v2) / (n - 1)
    plt.figure(figsize=(8, 5))
    
    v = np.linspace(min(v1, v2), max(v1, v2), 30)
    c = p1 * v1**n
    p = c / v**n
    
    plt.plot(v, p, 'k-', linewidth=3)
    plt.fill_between(v, p, color='red', alpha=0.3)
    plt.xlabel('Volume (m³)')
    plt.ylabel('Pressure (kPa)')
    plt.title(f'Polytropic Process (n={n})')
    plt.grid(True, alpha=0.3)
    plt.show()
    return round(work, 3)

# Sensible Heat Transfer
def q_sh(m, c, T1, T2):
    """
    Function for sensible heat transfer
    Input: mass (m), sp. heat (c), initial temp (T1), final temp (T2)
    Output: heat transfer
    """
    return m * c * (T2 - T1)

# Heat transfer during phase change
def q_phase_change(m, Ti, Tf, T_pc, c_bpc, c_apc, L):
    """
    Function for heat transfer during phase change
    Input: mass (m), initial temp (Ti), final temp (Tf),
    phase change temp (T_pc), sp. heat below/above phase change, latent heat (L)
    Output: heat interaction
    """
    if Ti > Tf:
        print('Process is cooling/freezing')
        return m * (c_bpc * (T_pc - Ti) - L + c_apc * (Tf - T_pc))
    else:
        print('Process is heating/melting')
        return m * (c_bpc * (T_pc - Ti) + L + c_apc * (Tf - T_pc))

print("Thermodynamic functions loaded successfully!")

Thermodynamic functions loaded successfully!

Example 1: Isothermal Compression

Air with mass 1.5 kg is compressed from 0.2 MPa to 1 MPa isothermally. Initial density is 1.21 kg/m³. Find work done

# Input Data
m = 1.50  # kg
p1 = 0.2E3  # kPa
p2 = 1.0E3  # kPa
density1 = 1.21  # kg/m^3

# Volume calculations
v1 = m / density1
v2 = (p1 * v1) / p2

# Calculate work for isothermal process
work = w_ct(p1, v1, v2)
print(f'Work output = {work} kJ')

Work output = -399.034 kJ

Example 2: Applying First Law

Using data from Example 1, if the system loses 5 kJ of heat, find the change in internal energy

# From previous example
work = -399.034  # kJ
heat = -5  # kJ (heat lost by system)

# First Law: Q = ?U + W
# Therefore: ?U = Q - W
delta_u = heat - work

print(f'Change in internal energy = {delta_u} kJ')
print(f'Internal energy increases by {abs(delta_u)} kJ')

Change in internal energy = 394.034 kJ
Internal energy increases by 394.034 kJ

Example 3: Phase Change Heat Transfer

Calculate heat removed when freezing 650 kg of tuna from 5°C to -12°C (freezing point: -2°C)

# Input data
m = 650  # kg
Ti = 5  # °C
Tf = -12  # °C  
T_pc = -2  # °C (freezing point)
c_bpc = 3.182  # kJ/kg-K (above freezing)
c_apc = 1.717  # kJ/kg-K (below freezing)
L = 234.5  # kJ/kg (latent heat of fusion)

# Calculate heat transfer
Q_net = q_phase_change(m, Ti, Tf, T_pc, c_bpc, c_apc, L)
print(f"Heat removed = {abs(Q_net)} kJ")

Process is cooling/freezing
Heat removed = 178063.6 kJ

Conclusion

Python provides powerful tools for modeling thermodynamic processes and applying the first law of thermodynamics. The functions demonstrated here can calculate work and heat interactions for various processes, enabling engineers to solve complex thermodynamic problems efficiently. The visualization capabilities help understand the relationship between pressure and volume in different processes.