Chapter 3: Key Catalytic Reactions

3.1 Haber-Bosch Process: Nitrogen Fixation

The Reaction

The Challenge

The N≡N triple bond is one of the strongest chemical bonds (945 kJ/mol). Breaking this bond requires:

• Iron catalyst with promoters (K₂O, Al₂O₃, CaO)
• High pressure (150-300 atm) to shift equilibrium toward products
• Moderate temperature (400-500°C) - compromise between kinetics and equilibrium

Code Example 1: Haber-Bosch Equilibrium Calculator

"""
Calculate Haber-Bosch equilibrium conversion at various conditions
Demonstrates the trade-off between temperature and pressure
"""
import numpy as np
import matplotlib.pyplot as plt

def haber_equilibrium(T_celsius, P_atm, initial_ratio=3):
    """
    Calculate equilibrium NH3 yield using simplified equilibrium expression
    N2 + 3H2 <=> 2NH3

    Uses van't Hoff equation approximation
    """
    T = T_celsius + 273.15  # Convert to Kelvin

    # Equilibrium constant (approximate, from thermodynamic data)
    # ln(Kp) ≈ 4710/T - 4.74 (simplified)
    ln_Kp = 4710 / T - 4.74
    Kp = np.exp(ln_Kp)

    # For stoichiometric feed (3H2:1N2), solving equilibrium
    # This is a simplified calculation
    # In reality, iterative solving is needed

    # Approximate conversion based on empirical correlation
    # Conversion increases with pressure, decreases with temperature
    x_eq = Kp * (P_atm / 100)**0.5 / (1 + Kp * (P_atm / 100)**0.5)
    x_eq = min(x_eq, 0.95)  # Cap at 95%

    return x_eq * 100  # Return as percentage

# Create meshgrid for pressure and temperature
temperatures = np.linspace(300, 600, 50)
pressures = [50, 100, 200, 300, 400]

# Plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Left: NH3 yield vs temperature at various pressures
colors = plt.cm.viridis(np.linspace(0, 1, len(pressures)))
for P, color in zip(pressures, colors):
    yields = [haber_equilibrium(T, P) for T in temperatures]
    ax1.plot(temperatures, yields, color=color, linewidth=2, label=f'{P} atm')

ax1.axvspan(400, 500, alpha=0.2, color='orange', label='Industrial range')
ax1.set_xlabel('Temperature (°C)', fontsize=12)
ax1.set_ylabel('NH₃ Equilibrium Yield (%)', fontsize=12)
ax1.set_title('Haber-Bosch: Temperature vs Yield', fontsize=13, fontweight='bold')
ax1.legend(title='Pressure', fontsize=9)
ax1.grid(alpha=0.3)
ax1.set_xlim(300, 600)
ax1.set_ylim(0, 100)

# Right: Reaction rate consideration
# Rate increases with T, but equilibrium decreases
rate_factor = np.exp(-(50000)/(8.314 * (temperatures + 273.15)))
rate_factor = rate_factor / rate_factor.max() * 100

equilibrium_yield = [haber_equilibrium(T, 200) for T in temperatures]
effective_rate = np.array(rate_factor) * np.array(equilibrium_yield) / 100

ax2.plot(temperatures, rate_factor, 'r--', linewidth=2, label='Reaction Rate')
ax2.plot(temperatures, equilibrium_yield, 'b--', linewidth=2, label='Equilibrium Yield')
ax2.plot(temperatures, effective_rate, 'g-', linewidth=3, label='Effective Production')

ax2.axvline(x=450, color='orange', linestyle=':', linewidth=2)
ax2.annotate('Optimal T ≈ 450°C', xy=(450, 50), fontsize=10, ha='left')

ax2.set_xlabel('Temperature (°C)', fontsize=12)
ax2.set_ylabel('Relative Value (%)', fontsize=12)
ax2.set_title('Kinetics vs Thermodynamics Trade-off', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("\nHaber-Bosch Process Impact:")
print("• Annual production: ~180 million tons NH₃")
print("• Feeds ~50% of world population through fertilizers")
print("• Consumes ~1.5% of global energy")
print("• Nobel Prize: Fritz Haber (1918), Carl Bosch (1931)")

3.2 Fischer-Tropsch Synthesis

Syngas to Liquid Fuels

Fischer-Tropsch (FT) synthesis converts synthesis gas (CO + H₂) into liquid hydrocarbons:

Anderson-Schulz-Flory Distribution

FT products follow a statistical distribution governed by chain growth probability α:

where $W_n$ is the weight fraction of hydrocarbon with n carbons.

Code Example 2: FT Product Distribution

"""
Calculate Anderson-Schulz-Flory distribution for Fischer-Tropsch synthesis
Shows how chain growth probability affects product distribution
"""
import numpy as np
import matplotlib.pyplot as plt

def asf_distribution(n, alpha):
    """Anderson-Schulz-Flory distribution"""
    return n * (1 - alpha)**2 * alpha**(n - 1)

# Carbon numbers
n = np.arange(1, 40)

# Different chain growth probabilities
alphas = [0.7, 0.8, 0.85, 0.9, 0.95]

# Create figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# Left: Weight distribution for different alpha
colors = plt.cm.plasma(np.linspace(0.2, 0.9, len(alphas)))
for alpha, color in zip(alphas, colors):
    W = asf_distribution(n, alpha)
    ax1.plot(n, W * 100, color=color, linewidth=2, label=f'α = {alpha}')

ax1.axvspan(5, 11, alpha=0.2, color='green', label='Gasoline (C5-C11)')
ax1.axvspan(12, 20, alpha=0.2, color='blue', label='Diesel (C12-C20)')
ax1.set_xlabel('Carbon Number', fontsize=12)
ax1.set_ylabel('Weight Fraction (%)', fontsize=12)
ax1.set_title('FT Product Distribution (ASF)', fontsize=13, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(alpha=0.3)
ax1.set_xlim(0, 40)

# Right: Selectivity to different fractions vs alpha
alphas_fine = np.linspace(0.5, 0.98, 100)
n_full = np.arange(1, 100)

gasoline_sel = []  # C5-C11
diesel_sel = []    # C12-C20
wax_sel = []       # C21+
methane_sel = []   # C1

for alpha in alphas_fine:
    W = asf_distribution(n_full, alpha)
    methane_sel.append(W[0] * 100)
    gasoline_sel.append(np.sum(W[4:11]) * 100)
    diesel_sel.append(np.sum(W[11:20]) * 100)
    wax_sel.append(np.sum(W[20:]) * 100)

ax2.plot(alphas_fine, methane_sel, 'r-', linewidth=2, label='Methane (C1)')
ax2.plot(alphas_fine, gasoline_sel, 'g-', linewidth=2, label='Gasoline (C5-C11)')
ax2.plot(alphas_fine, diesel_sel, 'b-', linewidth=2, label='Diesel (C12-C20)')
ax2.plot(alphas_fine, wax_sel, 'm-', linewidth=2, label='Wax (C21+)')

ax2.axvline(x=0.85, color='gray', linestyle='--', alpha=0.5)
ax2.annotate('Typical Co\nα ≈ 0.85', xy=(0.85, 40), fontsize=9, ha='center')

ax2.set_xlabel('Chain Growth Probability (α)', fontsize=12)
ax2.set_ylabel('Selectivity (%)', fontsize=12)
ax2.set_title('FT Selectivity vs Chain Growth', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9)
ax2.grid(alpha=0.3)

plt.tight_layout()
plt.show()

print("\nFischer-Tropsch Applications:")
print("• GTL (Gas-to-Liquids): Qatar, Malaysia, South Africa")
print("• CTL (Coal-to-Liquids): China, South Africa (Sasol)")
print("• BTL (Biomass-to-Liquids): Emerging renewable pathway")
print("• Typical conditions: 200-300°C, 10-40 bar, Co or Fe catalyst")

3.3 Catalytic Cracking and Reforming

Fluid Catalytic Cracking (FCC)

FCC is the largest volume catalytic process in the world, converting heavy oil fractions into gasoline:

• Catalyst: Zeolite Y + matrix (alumina, clay)
• Conditions: 500-550°C, ~1-3 bar
• Capacity: ~14 million barrels/day worldwide

Catalytic Reforming

Reforming converts low-octane naphtha into high-octane gasoline and aromatics:

• Catalyst: Pt/Re on acidic Al₂O₃
• Reactions: Dehydrogenation, isomerization, cyclization
• By-product: Hydrogen (valuable for hydroprocessing)

Code Example 3: Octane Number Enhancement

"""
Visualize catalytic reforming's effect on octane number
Shows transformation of feed to reformate
"""
import numpy as np
import matplotlib.pyplot as plt

# Hydrocarbon composition and octane numbers
# Before and after reforming

feed_components = {
    'n-Hexane': (40, 25),      # (wt%, RON)
    'n-Heptane': (25, 0),
    'Methylcyclopentane': (15, 91),
    'Cyclohexane': (10, 83),
    'Benzene': (5, 106),
    'Toluene': (5, 120),
}

reformate_components = {
    'n-Hexane': (5, 25),
    'n-Heptane': (3, 0),
    'Methylcyclopentane': (2, 91),
    'Cyclohexane': (5, 83),
    'Benzene': (25, 106),
    'Toluene': (35, 120),
    'Xylenes': (15, 117),
    'C9+ Aromatics': (10, 115),
}

def calculate_ron(components):
    """Calculate Research Octane Number from components"""
    total_wt = sum(v[0] for v in components.values())
    ron = sum(v[0] * v[1] for v in components.values()) / total_wt
    return ron

feed_ron = calculate_ron(feed_components)
reformate_ron = calculate_ron(reformate_components)

# Visualization
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(15, 5))

# Feed composition
feed_names = list(feed_components.keys())
feed_values = [v[0] for v in feed_components.values()]
colors_feed = ['#ffcccc' if v[1] < 50 else '#ccffcc' if v[1] > 100 else '#ffffcc'
               for v in feed_components.values()]
ax1.barh(feed_names, feed_values, color=colors_feed, edgecolor='black')
ax1.set_xlabel('Weight %', fontsize=11)
ax1.set_title(f'Naphtha Feed\nRON = {feed_ron:.0f}', fontsize=12, fontweight='bold')
ax1.set_xlim(0, 45)

# Reformate composition
ref_names = list(reformate_components.keys())
ref_values = [v[0] for v in reformate_components.values()]
colors_ref = ['#ffcccc' if v[1] < 50 else '#ccffcc' if v[1] > 100 else '#ffffcc'
              for v in reformate_components.values()]
ax2.barh(ref_names, ref_values, color=colors_ref, edgecolor='black')
ax2.set_xlabel('Weight %', fontsize=11)
ax2.set_title(f'Reformate Product\nRON = {reformate_ron:.0f}', fontsize=12, fontweight='bold')
ax2.set_xlim(0, 45)

# RON comparison
categories = ['Feed\nNaphtha', 'Reformate', 'Target\n(Premium)']
rons = [feed_ron, reformate_ron, 95]
colors_ron = ['#ff9999', '#99ff99', '#9999ff']
ax3.bar(categories, rons, color=colors_ron, edgecolor='black', width=0.6)
ax3.set_ylabel('Research Octane Number', fontsize=11)
ax3.set_title('Octane Number Enhancement', fontsize=12, fontweight='bold')
ax3.set_ylim(0, 130)
ax3.axhline(y=87, color='orange', linestyle='--', label='Regular (87)')
ax3.axhline(y=95, color='green', linestyle='--', label='Premium (95)')
ax3.legend(fontsize=9)

for i, v in enumerate(rons):
    ax3.text(i, v + 3, f'{v:.0f}', ha='center', fontsize=11, fontweight='bold')

plt.tight_layout()
plt.show()

print("\nReforming Process Summary:")
print(f"• Feed RON: {feed_ron:.0f} → Reformate RON: {reformate_ron:.0f}")
print(f"• Octane boost: +{reformate_ron - feed_ron:.0f} points")
print("• Key reactions: Dehydrocyclization, Isomerization, Aromatization")
print("• Hydrogen yield: 2-4% (valuable by-product)")

3.4 Hydrogenation Reactions

Types of Hydrogenation

Code Example 4: Hydrogenation Kinetics

"""
Model hydrogenation reaction kinetics
Compare different catalyst activities
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

def hydrogenation_kinetics(C, t, k, K_H, K_S, P_H2):
    """
    Langmuir-Hinshelwood kinetics for hydrogenation
    Rate = k * theta_H * theta_S
    where theta is surface coverage
    """
    C_substrate = C[0]

    # Surface coverages
    theta_H = K_H * P_H2 / (1 + K_H * P_H2 + K_S * C_substrate)
    theta_S = K_S * C_substrate / (1 + K_H * P_H2 + K_S * C_substrate)

    # Reaction rate
    rate = k * theta_H * theta_S

    return [-rate]  # Substrate consumption

# Time span
t = np.linspace(0, 60, 200)  # 60 minutes

# Initial substrate concentration
C0 = [1.0]  # mol/L

# Different catalysts with different parameters
catalysts = {
    'Pd/C': {'k': 0.5, 'K_H': 10, 'K_S': 5},
    'Pt/C': {'k': 0.3, 'K_H': 8, 'K_S': 8},
    'Ni/SiO₂': {'k': 0.1, 'K_H': 15, 'K_S': 3},
}

P_H2 = 1.0  # atm H2 pressure

# Solve and plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

colors = ['#e91e63', '#2196f3', '#4caf50']
for (name, params), color in zip(catalysts.items(), colors):
    solution = odeint(hydrogenation_kinetics, C0, t,
                      args=(params['k'], params['K_H'], params['K_S'], P_H2))

    conversion = (1 - solution[:, 0] / C0[0]) * 100
    ax1.plot(t, solution[:, 0], color=color, linewidth=2, label=name)
    ax2.plot(t, conversion, color=color, linewidth=2, label=name)

ax1.set_xlabel('Time (min)', fontsize=12)
ax1.set_ylabel('Substrate Concentration (mol/L)', fontsize=12)
ax1.set_title('Hydrogenation: Concentration Profile', fontsize=13, fontweight='bold')
ax1.legend(fontsize=10)
ax1.grid(alpha=0.3)

ax2.set_xlabel('Time (min)', fontsize=12)
ax2.set_ylabel('Conversion (%)', fontsize=12)
ax2.set_title('Hydrogenation: Conversion Profile', fontsize=13, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(alpha=0.3)
ax2.set_ylim(0, 100)

# Add t50 markers
ax2.axhline(y=50, color='gray', linestyle='--', alpha=0.5)
ax2.text(55, 52, 't₅₀', fontsize=10, color='gray')

plt.tight_layout()
plt.show()

# Calculate t50 for each catalyst
print("\nTime to 50% Conversion (t₅₀):")
for name, params in catalysts.items():
    solution = odeint(hydrogenation_kinetics, C0, t,
                      args=(params['k'], params['K_H'], params['K_S'], P_H2))
    conversion = (1 - solution[:, 0] / C0[0]) * 100
    t50_idx = np.argmin(np.abs(conversion - 50))
    print(f"  {name}: {t[t50_idx]:.1f} min")

3.5 Selective Oxidation

Partial vs Complete Oxidation

Selective oxidation produces valuable chemicals while avoiding complete combustion to CO₂:

3.6 Cross-Coupling Reactions (2010 Nobel Prize)

Major Cross-Coupling Reactions

Suzuki Coupling Mechanism

Code Example 5: Cross-Coupling Yield Optimization

"""
Analyze Suzuki coupling reaction conditions
Demonstrates effect of catalyst loading and temperature
"""
import numpy as np
import matplotlib.pyplot as plt

def suzuki_yield(temp_C, cat_loading_mol_pct, time_h, base_strength=1.0):
    """
    Empirical model for Suzuki coupling yield
    Based on typical reaction behavior
    """
    # Temperature effect (optimal around 80°C)
    temp_factor = np.exp(-((temp_C - 80) / 30)**2)

    # Catalyst loading effect (diminishing returns above 2%)
    cat_factor = 1 - np.exp(-cat_loading_mol_pct / 1.5)

    # Time effect (approaches completion)
    time_factor = 1 - np.exp(-time_h / 4)

    # Base yield around 90% under optimal conditions
    yield_pct = 95 * temp_factor * cat_factor * time_factor * base_strength

    return min(yield_pct, 99)

# Create analysis
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 1. Temperature effect
temps = np.linspace(20, 120, 100)
yields_temp = [suzuki_yield(T, 2.0, 12) for T in temps]
axes[0].plot(temps, yields_temp, 'b-', linewidth=2)
axes[0].axvline(x=80, color='green', linestyle='--', alpha=0.7, label='Optimal')
axes[0].set_xlabel('Temperature (°C)', fontsize=11)
axes[0].set_ylabel('Yield (%)', fontsize=11)
axes[0].set_title('Temperature Effect\n(2 mol% Pd, 12h)', fontsize=12, fontweight='bold')
axes[0].legend()
axes[0].grid(alpha=0.3)
axes[0].set_ylim(0, 100)

# 2. Catalyst loading effect
loadings = np.linspace(0.1, 5, 100)
yields_cat = [suzuki_yield(80, L, 12) for L in loadings]
axes[1].plot(loadings, yields_cat, 'r-', linewidth=2)
axes[1].axvline(x=2, color='green', linestyle='--', alpha=0.7, label='Standard')
axes[1].set_xlabel('Catalyst Loading (mol%)', fontsize=11)
axes[1].set_ylabel('Yield (%)', fontsize=11)
axes[1].set_title('Catalyst Loading Effect\n(80°C, 12h)', fontsize=12, fontweight='bold')
axes[1].legend()
axes[1].grid(alpha=0.3)
axes[1].set_ylim(0, 100)

# 3. Time course
times = np.linspace(0, 24, 100)
yields_time = [suzuki_yield(80, 2.0, t) for t in times]
axes[2].plot(times, yields_time, 'g-', linewidth=2)
axes[2].axhline(y=90, color='orange', linestyle='--', alpha=0.7, label='90% target')
axes[2].set_xlabel('Time (hours)', fontsize=11)
axes[2].set_ylabel('Yield (%)', fontsize=11)
axes[2].set_title('Reaction Progress\n(80°C, 2 mol% Pd)', fontsize=12, fontweight='bold')
axes[2].legend()
axes[2].grid(alpha=0.3)
axes[2].set_ylim(0, 100)

plt.tight_layout()
plt.show()

print("\nOptimal Suzuki Coupling Conditions:")
print("• Temperature: 60-100°C (commonly 80°C)")
print("• Catalyst: 1-5 mol% Pd (commonly 2 mol%)")
print("• Base: K₂CO₃, Cs₂CO₃, or KOH")
print("• Solvent: DMF, DME, THF/H₂O, or toluene/H₂O")
print("• Time: 6-24 hours for completion")

3.7 Polymerization Catalysis

Ziegler-Natta Catalysts (1963 Nobel Prize)

Ziegler-Natta catalysts enabled the production of stereoregular polymers, revolutionizing the plastics industry:

• Discovery: Karl Ziegler (polyethylene), Giulio Natta (polypropylene)
• Catalyst: TiCl₄ + AlR₃ (or supported TiCl₄/MgCl₂)
• Products: HDPE, isotactic PP, LLDPE
• Production: >150 million tons/year

Metallocene Catalysts

Metallocene catalysts are well-defined single-site catalysts that produce polymers with precise control:

• Structure: Cp₂MX₂ (Cp = cyclopentadienyl, M = Ti, Zr, Hf)
• Advantages: Narrow molecular weight distribution, controlled tacticity
• Activator: MAO (methylaluminoxane) or borates

3.8 Three-Way Catalytic Converter

Automotive Emission Control

The three-way catalyst (TWC) simultaneously converts three pollutants:

1. CO oxidation: 2CO + O₂ → 2CO₂
2. HC oxidation: CₓHᵧ + O₂ → CO₂ + H₂O
3. NOₓ reduction: 2NOₓ + 2CO → N₂ + 2CO₂

Code Example 6: TWC Light-off Behavior

"""
Model three-way catalyst light-off behavior
Shows conversion vs temperature for CO, HC, and NOx
"""
import numpy as np
import matplotlib.pyplot as plt

def twc_conversion(T, T50, steepness=0.1):
    """Sigmoid conversion profile for TWC"""
    return 100 / (1 + np.exp(-steepness * (T - T50)))

# Temperature range
T = np.linspace(100, 500, 200)

# Light-off temperatures (T50) for each pollutant
T50_CO = 200
T50_HC = 250
T50_NOx = 280

# Calculate conversions
conv_CO = twc_conversion(T, T50_CO, 0.08)
conv_HC = twc_conversion(T, T50_HC, 0.06)
conv_NOx = twc_conversion(T, T50_NOx, 0.05)

# Create plot
fig, ax = plt.subplots(figsize=(10, 6))

ax.plot(T, conv_CO, 'b-', linewidth=2.5, label='CO')
ax.plot(T, conv_HC, 'g-', linewidth=2.5, label='HC')
ax.plot(T, conv_NOx, 'r-', linewidth=2.5, label='NOₓ')

ax.axhline(y=90, color='gray', linestyle='--', alpha=0.5, label='90% target')
ax.axvspan(150, 300, alpha=0.1, color='orange', label='Light-off region')

ax.set_xlabel('Temperature (°C)', fontsize=12)
ax.set_ylabel('Conversion (%)', fontsize=12)
ax.set_title('Three-Way Catalyst Light-off Curves', fontsize=14, fontweight='bold')
ax.legend(fontsize=10, loc='lower right')
ax.grid(alpha=0.3)
ax.set_xlim(100, 500)
ax.set_ylim(0, 105)

# Add T50 annotations
ax.annotate(f'T₅₀(CO)={T50_CO}°C', xy=(T50_CO, 50), xytext=(T50_CO-50, 30),
            fontsize=9, arrowprops=dict(arrowstyle='->', color='blue'))
ax.annotate(f'T₅₀(HC)={T50_HC}°C', xy=(T50_HC, 50), xytext=(T50_HC+30, 30),
            fontsize=9, arrowprops=dict(arrowstyle='->', color='green'))

plt.tight_layout()
plt.show()

print("\nThree-Way Catalyst Specifications:")
print("• Catalyst: Pt/Pd/Rh on CeO₂-ZrO₂/Al₂O₃")
print("• Precious metal loading: 30-100 g/ft³")
print("• Operating temperature: 300-900°C")
print("• Conversion efficiency: >95% after light-off")
print("• Air/fuel ratio: λ = 1 (stoichiometric)")

3.9 Chapter Summary

Exercises

Exercise 1: Haber-Bosch Optimization

Calculate the equilibrium NH₃ yield at 450°C and 200 atm. What happens if you increase pressure to 400 atm?

Exercise 2: FT Product Analysis

For α = 0.85, calculate the weight fractions of C1-C4, C5-C11, C12-C20, and C21+ products.

Exercise 3: Reaction Design

Design a Suzuki coupling reaction to synthesize biphenyl from bromobenzene and phenylboronic acid. Specify catalyst, base, solvent, and conditions.

Exercise 4: TWC Analysis

Why must the air/fuel ratio be precisely controlled (λ ≈ 1) for a three-way catalyst to work efficiently?

Next Chapter

In Chapter 4: Characterization and Analysis, we will learn how to analyze catalyst properties using BET surface area, chemisorption, spectroscopy (XPS, FTIR), microscopy (TEM, SEM), and operando techniques.