Modelling Stirling and Ericsson Cycles in Python

The Stirling cycle and Ericsson cycle are important thermodynamic cycles used in heat engines. Python provides excellent tools for modeling these cycles using matplotlib and pandas to visualize the pressure-volume relationships and calculate state properties.

Stirling Cycle

The Stirling cycle consists of four processes: two reversible isochoric (constant volume) and two reversible isothermal (constant temperature) processes. The ideal regenerative Stirling cycle has the same efficiency as the Carnot cycle in the same temperature range.

Thermodynamic Calculations

The cycle involves these key calculations for each process:

• Process 1-2 (Isothermal compression): Volume decreases at constant temperature
• Process 2-3 (Isochoric heating): Pressure increases at constant volume
• Process 3-4 (Isothermal expansion): Volume increases at constant temperature
• Process 4-1 (Isochoric cooling): Pressure decreases at constant volume

Python Implementation

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

def stirling_cycle(p_min, p_max, v_max, r, gamma):
    """
    Model the Stirling cycle
    
    Parameters:
    p_min: minimum pressure (Pa)
    p_max: maximum pressure (Pa) 
    v_max: maximum volume (m³)
    r: compression ratio
    gamma: adiabatic exponent
    """
    
    # Process 1-2 (Isothermal compression)
    p1 = p_min
    v1 = v_max
    v2 = v1 / r
    c1 = p1 * v1  # Isothermal constant
    p2 = c1 / v2
    
    # Process 2-3 (Isochoric heating)
    p3 = p_max
    v3 = v2
    
    # Process 3-4 (Isothermal expansion)
    c2 = p3 * v3  # Isothermal constant
    v4 = v1
    p4 = c2 / v4
    
    # Create the plot
    plt.figure(figsize=(10, 8))
    
    # Process 1-2
    v_12 = np.linspace(v2, v1, 100)
    p_12 = c1 / v_12
    plt.plot(v_12, p_12/1000, 'r-', linewidth=2, label='1-2 Isothermal')
    
    # Process 2-3
    v_23 = np.full(100, v3)
    p_23 = np.linspace(p2, p3, 100)
    plt.plot(v_23, p_23/1000, 'b-', linewidth=2, label='2-3 Isochoric')
    
    # Process 3-4
    v_34 = np.linspace(v3, v4, 100)
    p_34 = c2 / v_34
    plt.plot(v_34, p_34/1000, 'g-', linewidth=2, label='3-4 Isothermal')
    
    # Process 4-1
    v_41 = np.full(100, v1)
    p_41 = np.linspace(p4, p1, 100)
    plt.plot(v_41, p_41/1000, 'm-', linewidth=2, label='4-1 Isochoric')
    
    # Add state points
    plt.plot([v1, v2, v3, v4], [p1/1000, p2/1000, p3/1000, p4/1000], 'ko', markersize=8)
    plt.text(v1, p1/1000-20, '1', fontsize=12, ha='center')
    plt.text(v2, p2/1000-20, '2', fontsize=12, ha='center')
    plt.text(v3, p3/1000+20, '3', fontsize=12, ha='center')
    plt.text(v4, p4/1000+20, '4', fontsize=12, ha='center')
    
    plt.title('Stirling Cycle P-V Diagram')
    plt.xlabel('Volume (m³)')
    plt.ylabel('Pressure (kPa)')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    # Return state data
    data = {
        'Pressure (Pa)': [p1, p2, p3, p4],
        'Volume (m³)': [v1, v2, v3, v4],
        'Process': ['Start', 'Compression', 'Heating', 'Expansion']
    }
    
    return pd.DataFrame(data, index=[1, 2, 3, 4])

# Example usage
result = stirling_cycle(p_min=1e5, p_max=20e5, v_max=0.5, r=5, gamma=1.4)
print(result)

   Pressure (Pa)  Volume (m³)     Process
1      100000.0         0.50       Start
2      500000.0         0.10  Compression
3     2000000.0         0.10     Heating
4      400000.0         0.50   Expansion

Ericsson Cycle

The Ericsson cycle consists of two reversible isobaric (constant pressure) and two reversible isothermal processes. Like the Stirling cycle, it has the same efficiency as the Carnot cycle when operating between the same temperature limits.

Python Implementation

def ericsson_cycle(p_min, p_max, v_max, r, gamma):
    """
    Model the Ericsson cycle
    
    Parameters:
    p_min: minimum pressure (Pa)
    p_max: maximum pressure (Pa)
    v_max: maximum volume (m³)
    r: compression ratio
    gamma: adiabatic exponent
    """
    
    # Process 1-2 (Isobaric compression)
    p1 = p_min
    p2 = p1
    v1 = v_max
    v2 = v1 / r
    
    # Process 2-3 (Isothermal heating)
    c1 = p2 * v2  # Isothermal constant
    p3 = p_max
    v3 = c1 / p3
    
    # Process 3-4 (Isobaric expansion)
    p4 = p3
    c2 = p1 * v1  # Isothermal constant for 4-1
    v4 = c2 / p4
    
    # Create the plot
    plt.figure(figsize=(10, 8))
    
    # Process 1-2
    v_12 = np.linspace(v2, v1, 100)
    p_12 = np.full(100, p2)
    plt.plot(v_12, p_12/1000, 'r-', linewidth=2, label='1-2 Isobaric')
    
    # Process 2-3
    v_23 = np.linspace(v3, v2, 100)
    p_23 = c1 / v_23
    plt.plot(v_23, p_23/1000, 'b-', linewidth=2, label='2-3 Isothermal')
    
    # Process 3-4
    v_34 = np.linspace(v3, v4, 100)
    p_34 = np.full(100, p4)
    plt.plot(v_34, p_34/1000, 'g-', linewidth=2, label='3-4 Isobaric')
    
    # Process 4-1
    v_41 = np.linspace(v4, v1, 100)
    p_41 = c2 / v_41
    plt.plot(v_41, p_41/1000, 'm-', linewidth=2, label='4-1 Isothermal')
    
    # Add state points
    plt.plot([v1, v2, v3, v4], [p1/1000, p2/1000, p3/1000, p4/1000], 'ko', markersize=8)
    plt.text(v1, p1/1000-50, '1', fontsize=12, ha='center')
    plt.text(v2, p2/1000-50, '2', fontsize=12, ha='center')
    plt.text(v3, p3/1000+50, '3', fontsize=12, ha='center')
    plt.text(v4, p4/1000+50, '4', fontsize=12, ha='center')
    
    plt.title('Ericsson Cycle P-V Diagram')
    plt.xlabel('Volume (m³)')
    plt.ylabel('Pressure (kPa)')
    plt.legend()
    plt.grid(True, alpha=0.3)
    plt.show()
    
    # Return state data
    data = {
        'Pressure (Pa)': [p1, p2, p3, p4],
        '