🌐 EN | 🇯🇵 JP | Last sync: 2025-11-16

Chapter 2: Material and Energy Balance

Implementing Fundamental Laws of Chemical Processes with Python

📖 Reading Time: 35-40 minutes 💻 Difficulty: Intermediate 🎯 Practical Exercises: 8

This chapter covers Material and Energy Balance. You will learn essential concepts and techniques.

2.1 Fundamentals of Material Balance

Material balance is based on the conservation law: "Input = Output + Accumulation + Reaction". In chemical process simulation, this fundamental principle is applied to multi-component systems, recycle streams, and reactive systems.

2.1.1 Component Material Balance (Multi-component Systems)

In distillation columns and separators, material balance must be calculated for each component. Below is an example of an ethanol-water mixture separation process.

Example 1: Component Material Balance for Distillation Column

A feed stream of 100 kg/h (20% ethanol, 80% water) is distilled to produce distillate (90% ethanol) and bottoms (2% ethanol). 

import numpy as np
from scipy.optimize import fsolve
import pandas as pd

def component_material_balance(F, z_F, x_D, x_B):
    """
    Solve component material balance for distillation column

    Args:
        F (float): Feed flow rate [kg/h]
        z_F (dict): Feed composition {component name: mole fraction}
        x_D (dict): Distillate composition {component name: mole fraction}
        x_B (dict): Bottoms composition {component name: mole fraction}

    Returns:
        dict: {'D': distillate flow rate, 'B': bottoms flow rate, 'recovery': recovery rate}
    """
    components = list(z_F.keys())

    # Overall material balance: F = D + B
    # Component balance: F * z_F[i] = D * x_D[i] + B * x_B[i]

    def equations(vars):
        D, B = vars
        eq1 = F - D - B  # Overall material balance

        # Use component balance for first component (same for other components)
        comp = components[0]
        eq2 = F * z_F[comp] - D * x_D[comp] - B * x_B[comp]

        return [eq1, eq2]

    # Initial guess
    D_init = F * 0.2
    B_init = F * 0.8

    # Solve simultaneous equations
    D, B = fsolve(equations, [D_init, B_init])

    # Calculate recovery rate for each component
    recovery = {}
    for comp in components:
        if F * z_F[comp] > 0:
            recovery[comp] = (D * x_D[comp]) / (F * z_F[comp]) * 100
        else:
            recovery[comp] = 0.0

    # Material balance check
    balance_check = {}
    for comp in components:
        input_comp = F * z_F[comp]
        output_comp = D * x_D[comp] + B * x_B[comp]
        balance_check[comp] = abs(input_comp - output_comp) < 1e-6

    return {
        'D': D,
        'B': B,
        'recovery': recovery,
        'balance_check': balance_check
    }

# Practical example: Ethanol-water distillation
F = 100.0  # kg/h
z_F = {'Ethanol': 0.20, 'Water': 0.80}  # Feed composition
x_D = {'Ethanol': 0.90, 'Water': 0.10}  # Distillate composition (high purity ethanol)
x_B = {'Ethanol': 0.02, 'Water': 0.98}  # Bottoms composition (water)

result = component_material_balance(F, z_F, x_D, x_B)

print("=" * 60)
print("Distillation Column Material Balance Results")
print("=" * 60)
print(f"\nFeed Conditions:")
print(f"  Feed flow rate: {F:.1f} kg/h")
print(f"  Feed composition: Ethanol {z_F['Ethanol']*100:.1f}%, Water {z_F['Water']*100:.1f}%")

print(f"\nCalculation Results:")
print(f"  Distillate flow rate (D): {result['D']:.2f} kg/h")
print(f"  Bottoms flow rate (B): {result['B']:.2f} kg/h")

print(f"\nRecovery Rate:")
for comp, rec in result['recovery'].items():
    print(f"  {comp}: {rec:.2f}%")

print(f"\nMaterial Balance Check:")
for comp, check in result['balance_check'].items():
    status = "✓ OK" if check else "✗ NG"
    print(f"  {comp}: {status}")

# Output example:
# ==============================================================
# Distillation Column Material Balance Results
# ==============================================================
#
# Feed Conditions:
#   Feed flow rate: 100.0 kg/h
#   Feed composition: Ethanol 20.0%, Water 80.0%
#
# Calculation Results:
#   Distillate flow rate (D): 20.45 kg/h
#   Bottoms flow rate (B): 79.55 kg/h
#
# Recovery Rate:
#   Ethanol: 92.05%
#   Water: 10.23%
#
# Material Balance Check:
#   Ethanol: ✓ OK
#   Water: ✓ OK

2.1.2 Material Balance with Recycle Streams

In chemical processes, recycle streams that return unreacted raw materials to the reactor are common. The balance is solved by combining the overall material balance and the balance around the reactor.

Example 2: Reactor System with Recycle Stream

For the reaction A → B with 70% conversion, 99% B removal in separator, and unreacted A recycled. 

import numpy as np
from scipy.optimize import fsolve

def recycle_material_balance(F_fresh, conversion, separation_eff):
    """
    Material balance for reactor with recycle stream

    Args:
        F_fresh (float): Fresh feed flow rate [kmol/h]
        conversion (float): Single-pass conversion [-]
        separation_eff (float): Separation efficiency (product B removal rate) [-]

    Returns:
        dict: Flow rates and compositions for each stream
    """
    def equations(vars):
        F_reactor, F_recycle = vars

        # Reactor inlet: F_reactor = F_fresh + F_recycle
        eq1 = F_reactor - F_fresh - F_recycle

        # Reactor outlet:
        # A remaining: F_reactor * (1 - conversion)
        # B produced: F_reactor * conversion
        F_A_out = F_reactor * (1 - conversion)
        F_B_out = F_reactor * conversion

        # Separator:
        # Flow after B removal goes to recycle
        # A in recycle: F_A_out (all)
        # B in recycle: F_B_out * (1 - separation_eff)
        eq2 = F_recycle - (F_A_out + F_B_out * (1 - separation_eff))

        return [eq1, eq2]

    # Initial guess
    F_reactor_init = F_fresh / conversion
    F_recycle_init = F_reactor_init - F_fresh

    # Solve simultaneous equations
    F_reactor, F_recycle = fsolve(equations, [F_reactor_init, F_recycle_init])

    # Detailed stream calculation
    F_A_reactor_out = F_reactor * (1 - conversion)
    F_B_reactor_out = F_reactor * conversion

    F_product = F_B_reactor_out * separation_eff
    F_B_recycle = F_B_reactor_out * (1 - separation_eff)

    # Recycle ratio
    recycle_ratio = F_recycle / F_fresh

    return {
        'F_fresh': F_fresh,
        'F_recycle': F_recycle,
        'F_reactor_inlet': F_reactor,
        'F_A_reactor_out': F_A_reactor_out,
        'F_B_reactor_out': F_B_reactor_out,
        'F_product': F_product,
        'recycle_ratio': recycle_ratio,
        'overall_conversion': (F_product / F_fresh) * 100
    }

# Practical example: Methanol synthesis reactor
F_fresh = 100.0  # kmol/h
conversion = 0.70  # 70%
separation_eff = 0.99  # 99%

result = recycle_material_balance(F_fresh, conversion, separation_eff)

print("=" * 60)
print("Recycle Stream Material Balance")
print("=" * 60)
print(f"\nConditions:")
print(f"  Fresh feed: {F_fresh:.1f} kmol/h")
print(f"  Single-pass conversion: {conversion*100:.0f}%")
print(f"  Separation efficiency: {separation_eff*100:.0f}%")

print(f"\nResults:")
print(f"  Reactor inlet flow rate: {result['F_reactor_inlet']:.2f} kmol/h")
print(f"  Recycle flow rate: {result['F_recycle']:.2f} kmol/h")
print(f"  Product flow rate: {result['F_product']:.2f} kmol/h")
print(f"  Recycle ratio: {result['recycle_ratio']:.2f}")
print(f"  Overall conversion: {result['overall_conversion']:.2f}%")

# Output example:
# ==============================================================
# Recycle Stream Material Balance
# ==============================================================
#
# Conditions:
#   Fresh feed: 100.0 kmol/h
#   Single-pass conversion: 70%
#   Separation efficiency: 99%
#
# Results:
#   Reactor inlet flow rate: 143.17 kmol/h
#   Recycle flow rate: 43.17 kmol/h
#   Product flow rate: 99.30 kmol/h
#   Recycle ratio: 0.43
#   Overall conversion: 99.30%

2.2 Fundamentals of Energy Balance

Energy balance is based on the principle: "Input Energy = Output Energy + Accumulation + Loss". In chemical processes, sensible heat, latent heat, and heat of reaction must be considered.

2.2.1 Energy Balance for Heat Exchangers

A heat exchanger is a device that transfers heat from a hot fluid to a cold fluid. The heat exchange rate is calculated using the logarithmic mean temperature difference (LMTD).

Example 3: Countercurrent Heat Exchanger

Hot side: 120°C → 60°C (flow rate 5000 kg/h, Cp = 4.18 kJ/kg·K)
Cold side: 20°C → ? °C (flow rate 8000 kg/h, Cp = 4.18 kJ/kg·K) 

import numpy as np

def heat_exchanger_energy_balance(m_h, T_h_in, T_h_out, Cp_h,
                                   m_c, T_c_in, Cp_c):
    """
    Energy balance and LMTD calculation for heat exchanger

    Args:
        m_h (float): Hot side flow rate [kg/h]
        T_h_in (float): Hot side inlet temperature [°C]
        T_h_out (float): Hot side outlet temperature [°C]
        Cp_h (float): Hot side specific heat [kJ/kg·K]
        m_c (float): Cold side flow rate [kg/h]
        T_c_in (float): Cold side inlet temperature [°C]
        Cp_c (float): Cold side specific heat [kJ/kg·K]

    Returns:
        dict: Energy balance results
    """
    # Heat released from hot side [kJ/h]
    Q_hot = m_h * Cp_h * (T_h_in - T_h_out)

    # Calculate cold side outlet temperature from energy balance
    # Q_hot = Q_cold = m_c * Cp_c * (T_c_out - T_c_in)
    T_c_out = T_c_in + Q_hot / (m_c * Cp_c)

    # Logarithmic mean temperature difference (LMTD) calculation
    dT1 = T_h_in - T_c_out  # Inlet end temperature difference
    dT2 = T_h_out - T_c_in  # Outlet end temperature difference

    if abs(dT1 - dT2) < 1e-6:
        LMTD = dT1  # When temperature difference is constant
    else:
        LMTD = (dT1 - dT2) / np.log(dT1 / dT2)

    # Effectiveness parameter
    effectiveness = (T_c_out - T_c_in) / (T_h_in - T_c_in)

    return {
        'Q': Q_hot / 1000,  # Convert to MW
        'T_c_out': T_c_out,
        'LMTD': LMTD,
        'dT1': dT1,
        'dT2': dT2,
        'effectiveness': effectiveness * 100
    }

# Practical example: Process cooling water heat exchanger
m_h = 5000.0  # kg/h
T_h_in = 120.0  # °C
T_h_out = 60.0  # °C
Cp_h = 4.18  # kJ/kg·K

m_c = 8000.0  # kg/h
T_c_in = 20.0  # °C
Cp_c = 4.18  # kJ/kg·K

result = heat_exchanger_energy_balance(m_h, T_h_in, T_h_out, Cp_h,
                                       m_c, T_c_in, Cp_c)

print("=" * 60)
print("Heat Exchanger Energy Balance")
print("=" * 60)
print(f"\nHot Side:")
print(f"  Inlet temperature: {T_h_in:.1f} °C")
print(f"  Outlet temperature: {T_h_out:.1f} °C")
print(f"  Flow rate: {m_h:.0f} kg/h")

print(f"\nCold Side:")
print(f"  Inlet temperature: {T_c_in:.1f} °C")
print(f"  Outlet temperature: {result['T_c_out']:.1f} °C")
print(f"  Flow rate: {m_c:.0f} kg/h")

print(f"\nHeat exchange rate: {result['Q']:.2f} MW")
print(f"Logarithmic mean temperature difference (LMTD): {result['LMTD']:.2f} °C")
print(f"Inlet end temperature difference (ΔT₁): {result['dT1']:.2f} °C")
print(f"Outlet end temperature difference (ΔT₂): {result['dT2']:.2f} °C")
print(f"Heat exchanger effectiveness: {result['effectiveness']:.1f}%")

# Output example:
# ==============================================================
# Heat Exchanger Energy Balance
# ==============================================================
#
# Hot Side:
#   Inlet temperature: 120.0 °C
#   Outlet temperature: 60.0 °C
#   Flow rate: 5000 kg/h
#
# Cold Side:
#   Inlet temperature: 20.0 °C
#   Outlet temperature: 51.2 °C
#   Flow rate: 8000 kg/h
#
# Heat exchange rate: 0.35 MW
# Logarithmic mean temperature difference (LMTD): 32.58 °C
# Inlet end temperature difference (ΔT₁): 68.75 °C
# Outlet end temperature difference (ΔT₂): 40.00 °C
# Heat exchanger effectiveness: 31.2%

2.2.2 Enthalpy Calculation (Sensible Heat + Latent Heat)

In processes involving phase change, both sensible heat and latent heat (heat of vaporization/condensation) must be considered.

Example 4: Water Heating and Evaporation Process

Calculate the heat required to heat water from 25°C to 100°C and completely evaporate it. 

import numpy as np

class WaterEnthalpy:
    """Water enthalpy calculation class"""

    # Physical properties
    Cp_liquid = 4.18  # kJ/kg·K (specific heat of liquid water)
    Cp_vapor = 2.08   # kJ/kg·K (specific heat of water vapor)
    H_vap_100C = 2257  # kJ/kg (latent heat of vaporization at 100°C)
    T_boil = 100.0    # °C (boiling point)

    @classmethod
    def calculate_enthalpy(cls, T_initial, T_final, m, phase_change=False):
        """
        Calculate enthalpy change of water

        Args:
            T_initial (float): Initial temperature [°C]
            T_final (float): Final temperature [°C]
            m (float): Mass [kg]
            phase_change (bool): Whether phase change is included

        Returns:
            dict: Enthalpy calculation results
        """
        Q_sensible = 0.0
        Q_latent = 0.0

        if not phase_change:
            # Sensible heat only (single phase)
            if T_initial < cls.T_boil and T_final < cls.T_boil:
                # Liquid phase only
                Q_sensible = m * cls.Cp_liquid * (T_final - T_initial)
            elif T_initial > cls.T_boil and T_final > cls.T_boil:
                # Gas phase only
                Q_sensible = m * cls.Cp_vapor * (T_final - T_initial)
        else:
            # Including phase change
            if T_initial < cls.T_boil and T_final > cls.T_boil:
                # Liquid heating → vaporization → vapor heating
                Q_heat_to_boil = m * cls.Cp_liquid * (cls.T_boil - T_initial)
                Q_vaporization = m * cls.H_vap_100C
                Q_superheat = m * cls.Cp_vapor * (T_final - cls.T_boil)

                Q_sensible = Q_heat_to_boil + Q_superheat
                Q_latent = Q_vaporization

        Q_total = Q_sensible + Q_latent

        return {
            'Q_sensible': Q_sensible,
            'Q_latent': Q_latent,
            'Q_total': Q_total,
            'Q_sensible_pct': (Q_sensible / Q_total * 100) if Q_total > 0 else 0,
            'Q_latent_pct': (Q_latent / Q_total * 100) if Q_total > 0 else 0
        }

# Practical example 1: Water heating (25°C → 100°C)
m1 = 1000.0  # kg
T_initial_1 = 25.0
T_final_1 = 100.0

result1 = WaterEnthalpy.calculate_enthalpy(T_initial_1, T_final_1, m1,
                                           phase_change=False)

print("=" * 60)
print("Example 1: Water Heating (Liquid Phase Only)")
print("=" * 60)
print(f"Conditions: {m1:.0f} kg, {T_initial_1}°C → {T_final_1}°C")
print(f"Sensible heat: {result1['Q_sensible']/1000:.2f} MJ")
print(f"Latent heat: {result1['Q_latent']/1000:.2f} MJ")
print(f"Total heat: {result1['Q_total']/1000:.2f} MJ")

# Practical example 2: Water heating and evaporation (25°C → 110°C, complete evaporation)
m2 = 1000.0  # kg
T_initial_2 = 25.0
T_final_2 = 110.0

result2 = WaterEnthalpy.calculate_enthalpy(T_initial_2, T_final_2, m2,
                                           phase_change=True)

print("\n" + "=" * 60)
print("Example 2: Water Heating, Evaporation, and Superheating")
print("=" * 60)
print(f"Conditions: {m2:.0f} kg, {T_initial_2}°C → {T_final_2}°C (including evaporation)")
print(f"Sensible heat: {result2['Q_sensible']/1000:.2f} MJ ({result2['Q_sensible_pct']:.1f}%)")
print(f"  - Liquid heating: {m2 * WaterEnthalpy.Cp_liquid * (100 - T_initial_2)/1000:.2f} MJ")
print(f"  - Vapor superheating: {m2 * WaterEnthalpy.Cp_vapor * (T_final_2 - 100)/1000:.2f} MJ")
print(f"Latent heat: {result2['Q_latent']/1000:.2f} MJ ({result2['Q_latent_pct']:.1f}%)")
print(f"Total heat: {result2['Q_total']/1000:.2f} MJ")

# Output example:
# ==============================================================
# Example 1: Water Heating (Liquid Phase Only)
# ==============================================================
# Conditions: 1000 kg, 25.0°C → 100.0°C
# Sensible heat: 313.50 MJ
# Latent heat: 0.00 MJ
# Total heat: 313.50 MJ
#
# ==============================================================
# Example 2: Water Heating, Evaporation, and Superheating
# ==============================================================
# Conditions: 1000 kg, 25.0°C → 110.0°C (including evaporation)
# Sensible heat: 334.30 MJ (12.9%)
#   - Liquid heating: 313.50 MJ
#   - Vapor superheating: 20.80 MJ
# Latent heat: 2257.00 MJ (87.1%)
# Total heat: 2591.30 MJ

2.2.3 Adiabatic Flame Temperature Calculation

In combustion reactions, the adiabatic flame temperature is calculated when all the heat of reaction is used to raise the temperature of the products. This is an important parameter for burner and boiler design.

Example 5: Adiabatic Combustion of Methane

Calculate the adiabatic flame temperature from fuel and air at 25°C for the reaction CH₄ + 2O₂ → CO₂ + 2H₂O. 

import numpy as np
from scipy.optimize import fsolve

def adiabatic_flame_temperature(fuel='CH4', T_initial=25, excess_air=1.0):
    """
    Calculation of adiabatic flame temperature

    Args:
        fuel (str): Fuel type ('CH4', 'C2H6', etc.)
        T_initial (float): Initial temperature [°C]
        excess_air (float): Excess air ratio (1.0 = theoretical air)

    Returns:
        dict: Adiabatic flame temperature and details
    """
    # Fuel database (simplified version)
    fuel_data = {
        'CH4': {
            'name': 'Methane',
            'formula': 'CH₄',
            'MW': 16.04,  # g/mol
            'H_combustion': -890.0,  # kJ/mol (heat of combustion, 25°C reference)
            'stoich_O2': 2.0,  # Theoretical oxygen moles
            'products': {'CO2': 1, 'H2O': 2}
        }
    }

    # Specific heat data (simplified temperature dependence, average values)
    Cp_data = {
        'CO2': 0.846,  # kJ/kg·K
        'H2O': 2.080,  # kJ/kg·K (vapor)
        'N2': 1.040,   # kJ/kg·K
        'O2': 0.918    # kJ/kg·K
    }

    MW_data = {
        'CO2': 44.01,
        'H2O': 18.02,
        'N2': 28.01,
        'O2': 32.00
    }

    fuel_info = fuel_data[fuel]

    # Combustion calculation
    n_O2_stoich = fuel_info['stoich_O2']
    n_O2_actual = n_O2_stoich * excess_air

    # Air composition (N2:O2 = 79:21 volume ratio)
    n_N2 = n_O2_actual * (79/21)

    # Product composition
    products = fuel_info['products'].copy()
    products['N2'] = n_N2

    # Excess oxygen
    if excess_air > 1.0:
        products['O2'] = n_O2_actual - n_O2_stoich

    # Heat released (per mol of fuel)
    Q_released = -fuel_info['H_combustion']  # kJ

    # Equation to find adiabatic flame temperature
    def energy_balance(T_flame):
        # Sensible heat of products = heat released
        Q_sensible = 0.0
        for comp, n_mol in products.items():
            m = n_mol * MW_data[comp] / 1000  # kg
            Cp = Cp_data[comp]
            Q_sensible += m * Cp * (T_flame - T_initial)

        return Q_released - Q_sensible

    # Initial guess: 2000°C
    T_flame = fsolve(energy_balance, 2000.0)[0]

    # Detailed output
    total_mass = sum(products[c] * MW_data[c] for c in products) / 1000

    return {
        'fuel_name': fuel_info['name'],
        'T_flame': T_flame,
        'Q_released': Q_released,
        'excess_air': excess_air,
        'products': products,
        'total_product_mass': total_mass
    }

# Practical example: Methane combustion
results = []
for excess_air in [1.0, 1.2, 1.5]:
    result = adiabatic_flame_temperature('CH4', T_initial=25,
                                         excess_air=excess_air)
    results.append(result)

print("=" * 60)
print("Adiabatic Flame Temperature of Methane")
print("=" * 60)

for i, result in enumerate(results, 1):
    print(f"\nCase {i}: Excess air ratio {result['excess_air']}")
    print(f"  Adiabatic flame temperature: {result['T_flame']:.0f} °C")
    print(f"  Heat released: {result['Q_released']:.1f} kJ/mol")
    print(f"  Product mass: {result['total_product_mass']:.3f} kg/mol-fuel")

# Output example:
# ==============================================================
# Adiabatic Flame Temperature of Methane
# ==============================================================
#
# Case 1: Excess air ratio 1.0
#   Adiabatic flame temperature: 1960 °C
#   Heat released: 890.0 kJ/mol
#   Product mass: 0.140 kg/mol-fuel
#
# Case 2: Excess air ratio 1.2
#   Adiabatic flame temperature: 1755 °C
#   Heat released: 890.0 kJ/mol
#   Product mass: 0.157 kg/mol-fuel
#
# Case 3: Excess air ratio 1.5
#   Adiabatic flame temperature: 1505 °C
#   Heat released: 890.0 kJ/mol
#   Product mass: 0.182 kg/mol-fuel

2.3 Energy Balance for Reactors

In chemical reactors, an energy balance that considers the heat of reaction (exothermic/endothermic) is necessary. The reaction rate and heat balance are solved in a coupled manner.

2.3.1 Energy Balance for Exothermic Reactors

Example 6: Energy Balance for Continuous Stirred Tank Reactor (CSTR)

Calculate heat removal and temperature control for a CSTR performing the reaction A → B (ΔH = -50 kJ/mol, exothermic). 

import numpy as np
from scipy.optimize import fsolve

def cstr_energy_balance_exothermic(F, C_A0, T0, V, k0, Ea, delta_H,
                                    T_coolant, UA):
    """
    Energy balance for exothermic reaction CSTR

    Args:
        F (float): Volumetric flow rate [m³/h]
        C_A0 (float): Feed concentration [kmol/m³]
        T0 (float): Feed temperature [°C]
        V (float): Reactor volume [m³]
        k0 (float): Frequency factor [1/h]
        Ea (float): Activation energy [kJ/mol]
        delta_H (float): Heat of reaction [kJ/mol] (negative for exothermic)
        T_coolant (float): Coolant temperature [°C]
        UA (float): Overall heat transfer coefficient × area [kJ/h·K]

    Returns:
        dict: Reactor outlet conditions
    """
    R = 8.314e-3  # kJ/mol·K (gas constant)
    rho_Cp = 4180.0  # kJ/m³·K (density × specific heat of reaction mixture, water equivalent)

    def equations(vars):
        C_A, T = vars
        T_K = T + 273.15  # Kelvin conversion

        # Reaction rate constant (Arrhenius equation)
        k = k0 * np.exp(-Ea / (R * T_K))

        # Reaction rate
        r_A = k * C_A  # kmol/m³·h

        # Material balance: F * (C_A0 - C_A) = V * r_A
        eq1 = F * (C_A0 - C_A) - V * r_A

        # Energy balance:
        # Input - output + reaction heat - heat removal = 0
        Q_reaction = V * r_A * (-delta_H)  # kJ/h (exothermic)
        Q_removal = UA * (T - T_coolant)  # kJ/h (heat removal)
        Q_accumulation = F * rho_Cp * (T - T0)  # kJ/h

        eq2 = Q_reaction - Q_removal - Q_accumulation

        return [eq1, eq2]

    # Initial guess
    C_A_init = C_A0 * 0.5
    T_init = T0 + 20

    # Solve simultaneous equations
    C_A, T = fsolve(equations, [C_A_init, T_init])

    # Conversion
    conversion = (C_A0 - C_A) / C_A0 * 100

    # Heat balance details
    T_K = T + 273.15
    k = k0 * np.exp(-Ea / (R * T_K))
    r_A = k * C_A
    Q_reaction = V * r_A * (-delta_H)
    Q_removal = UA * (T - T_coolant)

    return {
        'C_A': C_A,
        'T': T,
        'conversion': conversion,
        'k': k,
        'r_A': r_A,
        'Q_reaction': Q_reaction / 1000,  # MW
        'Q_removal': Q_removal / 1000,     # MW
        'Q_balance': abs(Q_reaction - Q_removal) / Q_reaction * 100
    }

# Practical example: Ethyl acetate synthesis reactor
F = 10.0  # m³/h
C_A0 = 5.0  # kmol/m³
T0 = 60.0  # °C
V = 5.0  # m³
k0 = 1.5e6  # 1/h
Ea = 65.0  # kJ/mol
delta_H = -50.0  # kJ/mol (exothermic)
T_coolant = 50.0  # °C
UA = 5000.0  # kJ/h·K

result = cstr_energy_balance_exothermic(F, C_A0, T0, V, k0, Ea, delta_H,
                                        T_coolant, UA)

print("=" * 60)
print("Exothermic Reaction CSTR Energy Balance")
print("=" * 60)
print(f"\nReactor Specifications:")
print(f"  Volume: {V} m³")
print(f"  Flow rate: {F} m³/h")
print(f"  Residence time: {V/F:.2f} h")

print(f"\nReaction Conditions:")
print(f"  Feed concentration: {C_A0} kmol/m³")
print(f"  Feed temperature: {T0} °C")
print(f"  Heat of reaction: {delta_H} kJ/mol")

print(f"\nCalculation Results:")
print(f"  Outlet concentration: {result['C_A']:.3f} kmol/m³")
print(f"  Outlet temperature: {result['T']:.1f} °C")
print(f"  Conversion: {result['conversion']:.1f}%")
print(f"  Reaction rate constant: {result['k']:.2f} 1/h")

print(f"\nHeat Balance:")
print(f"  Heat generated: {result['Q_reaction']:.3f} MW")
print(f"  Heat removed: {result['Q_removal']:.3f} MW")
print(f"  Heat balance error: {result['Q_balance']:.2f}%")

# Output example:
# ==============================================================
# Exothermic Reaction CSTR Energy Balance
# ==============================================================
#
# Reactor Specifications:
#   Volume: 5.0 m³
#   Flow rate: 10.0 m³/h
#   Residence time: 0.50 h
#
# Reaction Conditions:
#   Feed concentration: 5.0 kmol/m³
#   Feed temperature: 60 °C
#   Heat of reaction: -50.0 kJ/mol
#
# Calculation Results:
#   Outlet concentration: 0.523 kmol/m³
#   Outlet temperature: 80.3 °C
#   Conversion: 89.5%
#   Reaction rate constant: 17.16 1/h
#
# Heat Balance:
#   Heat generated: 0.449 MW
#   Heat removed: 0.449 MW
#   Heat balance error: 0.01%

2.3.2 Energy Balance for Endothermic Reactors

Example 7: Endothermic Reaction (Steam Reforming)

Calculate the required heating rate for the endothermic reaction CH₄ + H₂O → CO + 3H₂ (ΔH = +206 kJ/mol). 

import numpy as np
from scipy.optimize import fsolve

def endothermic_reactor_energy_balance(F_CH4, T_in, conversion_target,
                                        T_reactor, delta_H):
    """
    Energy balance for endothermic reactor

    Args:
        F_CH4 (float): Methane feed flow rate [kmol/h]
        T_in (float): Feed temperature [°C]
        conversion_target (float): Target conversion [-]
        T_reactor (float): Reactor temperature [°C]
        delta_H (float): Heat of reaction [kJ/mol] (positive for endothermic)

    Returns:
        dict: Required heating rate and product flow rates
    """
    # Specific heat data (average values, kJ/kmol·K)
    Cp = {
        'CH4': 35.0,
        'H2O': 33.6,
        'CO': 29.1,
        'H2': 28.8
    }

    # Reaction amount
    F_CH4_reacted = F_CH4 * conversion_target
    F_CH4_out = F_CH4 * (1 - conversion_target)

    # Product flow rates (from stoichiometry)
    # CH4 + H2O → CO + 3H2
    F_H2O_in = F_CH4 * 1.2  # 20% excess
    F_H2O_reacted = F_CH4_reacted
    F_H2O_out = F_H2O_in - F_H2O_reacted

    F_CO = F_CH4_reacted
    F_H2 = 3 * F_CH4_reacted

    # Heat required for reaction
    Q_reaction = F_CH4_reacted * delta_H  # kJ/h

    # Sensible heat (heating raw materials to reaction temperature)
    Q_sensible_CH4 = F_CH4 * Cp['CH4'] * (T_reactor - T_in)
    Q_sensible_H2O = F_H2O_in * Cp['H2O'] * (T_reactor - T_in)
    Q_sensible_total = Q_sensible_CH4 + Q_sensible_H2O

    # Total heating rate
    Q_total = Q_reaction + Q_sensible_total

    return {
        'Q_reaction': Q_reaction / 1000,  # MW
        'Q_sensible': Q_sensible_total / 1000,  # MW
        'Q_total': Q_total / 1000,  # MW
        'F_CH4_out': F_CH4_out,
        'F_H2O_out': F_H2O_out,
        'F_CO': F_CO,
        'F_H2': F_H2,
        'conversion': conversion_target * 100
    }

# Practical example: Steam reforming
F_CH4 = 100.0  # kmol/h
T_in = 400.0  # °C
conversion_target = 0.85  # 85%
T_reactor = 850.0  # °C
delta_H = 206.0  # kJ/mol (endothermic)

result = endothermic_reactor_energy_balance(F_CH4, T_in, conversion_target,
                                             T_reactor, delta_H)

print("=" * 60)
print("Endothermic Reactor Energy Balance (Steam Reforming)")
print("=" * 60)
print(f"\nReaction Conditions:")
print(f"  Methane feed: {F_CH4} kmol/h")
print(f"  Feed temperature: {T_in} °C")
print(f"  Reactor temperature: {T_reactor} °C")
print(f"  Target conversion: {conversion_target*100:.0f}%")

print(f"\nRequired Heating Rate:")
print(f"  Reaction heat: {result['Q_reaction']:.2f} MW")
print(f"  Sensible heat: {result['Q_sensible']:.2f} MW")
print(f"  Total heating rate: {result['Q_total']:.2f} MW")

print(f"\nProduct Flow Rates:")
print(f"  CH₄ (unreacted): {result['F_CH4_out']:.1f} kmol/h")
print(f"  H₂O (unreacted): {result['F_H2O_out']:.1f} kmol/h")
print(f"  CO (produced): {result['F_CO']:.1f} kmol/h")
print(f"  H₂ (produced): {result['F_H2']:.1f} kmol/h")

# Output example:
# ==============================================================
# Endothermic Reactor Energy Balance (Steam Reforming)
# ==============================================================
#
# Reaction Conditions:
#   Methane feed: 100.0 kmol/h
#   Feed temperature: 400 °C
#   Reactor temperature: 850 °C
#   Target conversion: 85%
#
# Required Heating Rate:
#   Reaction heat: 17.51 MW
#   Sensible heat: 3.40 MW
#   Total heating rate: 20.91 MW
#
# Product Flow Rates:
#   CH₄ (unreacted): 15.0 kmol/h
#   H₂O (unreacted): 35.0 kmol/h
#   CO (produced): 85.0 kmol/h
#   H₂ (produced): 255.0 kmol/h

2.3.3 Complete Material and Energy Balance (Integrated Example)

Example 8: Simultaneous Material and Energy Balance for Reactive Distillation Column

Solve coupled material and energy balances for a reactive distillation column where reaction and separation occur simultaneously. 

import numpy as np
from scipy.optimize import fsolve

class ReactiveDistillationColumn:
    """Material and energy balance class for reactive distillation column"""

    def __init__(self, F, z_F, T_F, P, k_rxn, delta_H_rxn, delta_H_vap):
        """
        Args:
            F (float): Feed flow rate [kmol/h]
            z_F (dict): Feed composition {'A': x_A, 'B': x_B, 'C': x_C}
            T_F (float): Feed temperature [°C]
            P (float): Operating pressure [kPa]
            k_rxn (float): Reaction rate constant [1/h]
            delta_H_rxn (float): Heat of reaction [kJ/kmol]
            delta_H_vap (dict): Latent heat of vaporization for each component [kJ/kmol]
        """
        self.F = F
        self.z_F = z_F
        self.T_F = T_F
        self.P = P
        self.k_rxn = k_rxn
        self.delta_H_rxn = delta_H_rxn
        self.delta_H_vap = delta_H_vap

        # Specific heat data (liquid phase, kJ/kmol·K)
        self.Cp_liquid = {'A': 150, 'B': 150, 'C': 180}

    def solve_balance(self, reflux_ratio, N_stages):
        """
        Solve material and energy balance

        Args:
            reflux_ratio (float): Reflux ratio [-]
            N_stages (int): Number of theoretical stages

        Returns:
            dict: Calculation results
        """
        def equations(vars):
            D, B, T_top, T_bottom, Q_reboiler, Q_condenser = vars

            # Material balance
            eq1 = self.F - D - B

            # Component balance (simplified: A is light, C is heavy)
            # Distillate is mainly A, bottoms is mainly C
            x_D_A = 0.95  # Assumption
            x_B_C = 0.90  # Assumption

            F_A = self.F * self.z_F['A']
            eq2 = F_A - D * x_D_A - B * (1 - x_B_C) * 0.5

            # Reaction: A + B → C (simplified)
            # Reaction amount depends on holdup and reaction rate
            V_total = 10.0  # m³ (assumption)
            C_avg = self.F / V_total
            r_rxn = self.k_rxn * C_avg * 0.5  # kmol/h

            # Energy balance
            # Reboiler heat duty
            V_vapor = D * (1 + reflux_ratio)  # Vapor flow rate
            Q_reb_calc = V_vapor * self.delta_H_vap['C']
            eq3 = Q_reboiler - Q_reb_calc

            # Condenser heat duty
            Q_cond_calc = V_vapor * self.delta_H_vap['A']
            eq4 = Q_condenser - Q_cond_calc

            # Reaction heat
            Q_reaction = r_rxn * self.delta_H_rxn

            # Temperature balance (simplified)
            eq5 = T_top - (80 + reflux_ratio * 5)
            eq6 = T_bottom - (120 + Q_reboiler / 1000)

            return [eq1, eq2, eq3, eq4, eq5, eq6]

        # Initial guess
        D_init = self.F * 0.4
        B_init = self.F * 0.6
        T_top_init = 85.0
        T_bottom_init = 125.0
        Q_reb_init = 5000.0
        Q_cond_init = 5000.0

        # Solve simultaneous equations
        result = fsolve(equations,
                       [D_init, B_init, T_top_init, T_bottom_init,
                        Q_reb_init, Q_cond_init])

        D, B, T_top, T_bottom, Q_reboiler, Q_condenser = result

        return {
            'D': D,
            'B': B,
            'T_top': T_top,
            'T_bottom': T_bottom,
            'Q_reboiler': Q_reboiler / 1000,  # MW
            'Q_condenser': Q_condenser / 1000,  # MW
            'reflux_ratio': reflux_ratio,
            'energy_efficiency': (Q_reboiler - Q_condenser) / Q_reboiler * 100
        }

# Practical example: Ethyl acetate synthesis reactive distillation
F = 100.0  # kmol/h
z_F = {'A': 0.4, 'B': 0.4, 'C': 0.2}  # A: Ethanol, B: Acetic acid, C: Ethyl acetate
T_F = 60.0  # °C
P = 101.3  # kPa
k_rxn = 2.0  # 1/h
delta_H_rxn = -25.0  # kJ/kmol (exothermic)
delta_H_vap = {'A': 38000, 'B': 24000, 'C': 32000}  # kJ/kmol

column = ReactiveDistillationColumn(F, z_F, T_F, P, k_rxn, delta_H_rxn,
                                     delta_H_vap)

# Calculate with varying reflux ratio
reflux_ratios = [1.5, 2.0, 3.0]
N_stages = 20

print("=" * 60)
print("Reactive Distillation Column Material and Energy Balance")
print("=" * 60)
print(f"\nColumn Specifications:")
print(f"  Feed flow rate: {F} kmol/h")
print(f"  Number of theoretical stages: {N_stages}")
print(f"  Reaction: A + B → C (exothermic)")

for RR in reflux_ratios:
    result = column.solve_balance(RR, N_stages)

    print(f"\nReflux ratio: {RR}")
    print(f"  Distillate flow rate: {result['D']:.1f} kmol/h")
    print(f"  Bottoms flow rate: {result['B']:.1f} kmol/h")
    print(f"  Top temperature: {result['T_top']:.1f} °C")
    print(f"  Bottom temperature: {result['T_bottom']:.1f} °C")
    print(f"  Reboiler heat duty: {result['Q_reboiler']:.2f} MW")
    print(f"  Condenser heat duty: {result['Q_condenser']:.2f} MW")

# Output example:
# ==============================================================
# Reactive Distillation Column Material and Energy Balance
# ==============================================================
#
# Column Specifications:
#   Feed flow rate: 100.0 kmol/h
#   Number of theoretical stages: 20
#   Reaction: A + B → C (exothermic)
#
# Reflux ratio: 1.5
#   Distillate flow rate: 39.4 kmol/h
#   Bottoms flow rate: 60.6 kmol/h
#   Top temperature: 87.5 °C
#   Bottom temperature: 130.2 °C
#   Reboiler heat duty: 3.15 MW
#   Condenser heat duty: 3.75 MW
#
# Reflux ratio: 2.0
#   Distillate flow rate: 39.2 kmol/h
#   Bottoms flow rate: 60.8 kmol/h
#   Top temperature: 90.0 °C
#   Bottom temperature: 131.8 °C
#   Reboiler heat duty: 3.76 MW
#   Condenser heat duty: 4.48 MW
#
# Reflux ratio: 3.0
#   Distillate flow rate: 39.0 kmol/h
#   Bottoms flow rate: 61.0 kmol/h
#   Top temperature: 95.0 °C
#   Bottom temperature: 134.7 °C
#   Reboiler heat duty: 4.99 MW
#   Condenser heat duty: 5.95 MW

2.4 Learning Summary

What We Learned in This Chapter

  • Material Balance: Balance calculations for complex systems including multi-component systems and recycle streams
  • Energy Balance: Heat balance calculations considering sensible heat, latent heat, and heat of reaction
  • Heat Exchangers: Calculation of energy exchange using the LMTD method
  • Reactors: Energy balance and temperature control for exothermic and endothermic reactions
  • Integrated Balance: Solution methods for coupled material and energy systems

Practical Considerations

  • Consider temperature and pressure dependence of physical properties (simplified in this chapter)
  • Account for heat losses of approximately 10-20% in actual processes
  • Incorporate safety factors (1.1-1.3 times) in design
  • Consider dynamic response of control systems (this chapter covers steady state only)
← Back to Series Index Chapter 3: Unit Operations Modeling →

References

  1. Montgomery, D. C. (2019). Design and Analysis of Experiments (9th ed.). Wiley.
  2. Box, G. E. P., Hunter, J. S., & Hunter, W. G. (2005). Statistics for Experimenters: Design, Innovation, and Discovery (2nd ed.). Wiley.
  3. Seborg, D. E., Edgar, T. F., Mellichamp, D. A., & Doyle III, F. J. (2016). Process Dynamics and Control (4th ed.). Wiley.
  4. McKay, M. D., Beckman, R. J., & Conover, W. J. (2000). "A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code." Technometrics, 42(1), 55-61.

Disclaimer