“Wave Tracker” for Pressure Swing Adsorption

Introduction

Pressure Swing Adsorption (PSA) is a widely adopted process design for gas separation. A foundational analytical framework for PSA is the equilibrium theory, first introduced by works such as Knaebel and Hill (1985). Equilibrium theory simplifies PSA analysis by neglecting mass transfer resistances, dispersive effects, and other dissipative phenomena. Under this framework, several assumptions are typically made:

• Linear adsorption isotherms.
• Instantaneous equilibrium between gas and adsorbed phases.
• Isothermal operation.
• Ideal gas behavior.

With these assumptions, the governing equations simplify significantly, facilitating quick numerical or even analytical solutions. Furthermore, the complex dynamics within adsorption columns can usually be adequately represented by a few concentration (or partial pressure) fronts and characteristic lines between these fronts, enhancing interpretability and visualization. Consequently, equilibrium theory provides a computationally inexpensive yet insightful approach to elucidate key characteristics of PSA processes across various boundary conditions, without resorting to computationally intensive PDE-based simulations. Equilibrium theory is particularly beneficial for:

• Guiding cycle designs and determining step sequences.
• Providing fast approximations, initial guesses, or surrogate models for more rigorous PDE-based analyses.
• Enabling simplified optimization and parametric sensitivity studies.
• Estimating process-level performance metrics for large-scale adsorbent screening and optimization.

In this blog post, I first provide an overview of the equilibrium theory as applied to PSA. Next, I demonstrate a Python code for wave tracking based on equilibrium theory. Lastly, I present several illustrative examples, showcasing the practical application of the code and interpreting the resulting wave diagrams.

Equilibrium Theory for PSA

Basic Setting

Consider a packed bed of length $L$ used for separating a binary gas mixture consisting of:

• Component $A$ (the heavy or more strongly adsorbed species),
• Component $B$ (the light or less strongly adsorbed species).

The total pressure is denoted by $P$, so $p_A+p_B=P$, where $p_A$ and $p_B$ are the partial pressures of $A$ and $B$, respectively.

Key assumptions:

1. Linear isotherms of the form: $$q_A=k_A{p}_A,q_B=k_B{p}_B,$$ where $q_A$ and $q_B$ are the adsorbed amounts of $A$ and $B$.
   
2. Instantaneous equilibrium between the gas and adsorbed phases (no mass-transfer resistance or kinetic limitations).
   
3. Ideal gas behavior and isothermal operation (temperature $T$ is constant).
   
4. The adsorbent bed has an interstitial void fraction $\epsilon$, meaning a fraction $\epsilon$ of the bed volume is accessible to gas flow, and $(1-\epsilon )$ is occupied by the solid adsorbent.
   
5. Under these assumptions, the evolution of gas composition and pressure along the bed can be described by a system of hyperbolic PDEs (or derived characteristic equations) which are simplified by the equilibrium and linear isotherm assumptions.

Dimensionless Groups

To simplify the governing equations, it is helpful to define: $${\beta }_A=\frac{1}{1+\frac{(1-\epsilon )}{\epsilon }{k}_A},{\beta }_B=\frac{1}{1+\frac{(1-\epsilon )}{\epsilon }{k}_B},\beta =\frac{{\beta }_B}{{\beta }_A}\le 1.$$ These dimensionless parameters incorporate the adsorption constants $k_A$, $k_B$ and the void fraction $\epsilon$. Notably, $\beta$ provides a measure of the separation factor: for large $k_A$ relative to $k_B$, $\beta$ becomes small, reflecting stronger selectivity for the heavy species $A$.

Shocks and Waves

In a typical PSA cycle, the column may experience phases of constant pressure (e.g., during a high-pressure feed step) and phases of changing pressure (e.g., blowdown or pressurization steps). During these steps, characteristic curves in the $(t,z)$ plane can converge (forming shocks) or diverge (forming expansion or simple waves). I outline the resulting wave phenomena below.

Constant-Pressure Operation

Suppose the gas enters the bed at one end at a fixed total pressure $P$. For example, during a high-pressure feed step, the velocity $u$ in the bed generally depends on the local composition $y$ (the mole fraction of the heavy component $A$ in the gas phase). A common result under constant $P$ is: $$u(z)\propto \frac{1}{1+(\beta -1)y(z)}.$$

• Shock waves form when a gas richer in the heavy component $A$ pushes into a region with lower $A$ concentration. In this case, characteristics merge, and the shock speed $u_S$ (in dimensionless form) is given by $$u_S={\beta }_A\frac{u_2{y}_2-u_1{y}_1}{y_2-y_1},$$ where $(u_1,y_1)$ and $(u_2,y_2)$ are the velocity and composition just ahead of and behind the shock, respectively.
• Simple waves (or expansion waves) appear when a more dilute gas (lower $y$) contacts a region of higher $y$. In this case, the characteristics diverge. Mathematically, the velocity and concentration along the characteristic can be determined by integrating the relevant PDE or by directly applying the method of characteristics, which yields a continuous fan of solutions rather than a jump (shock).

Changing-Pressure Steps

In other steps (e.g., blowdown or pressurization), the total pressure $P$ varies with time. One end of the column is typically closed (no flow), while gas either exits or enters from the other end. In this scenario, the local velocity $u$ and the bed composition are governed by relationships such as $$\frac{1}{{\beta }_B}\frac{\partial P}{\partial t}+\dots =0,\text{and}u=\text{function}(P,\frac{dP}{dt},z),$$ derived from the continuity equations and the assumption of instantaneous equilibrium. By integrating along characteristic curves, one obtains explicit expressions for how composition $y$ changes with $P$, and how wave fronts (shocks or expansions) propagate through the column. The shock speed in a variable-pressure step is similarly given by $$u_S={\beta }_A\frac{u_2{y}_2-u_1{y}_1}{y_2-y_1},$$ but with $u_1$ and $u_2$ determined by the pressure-change relationships (often involving $\frac{dP}{dt}$ and the dimensionless groups ${\beta }_A,{\beta }_B$).

Under these linear and equilibrium assumptions, the bed dynamics during a PSA cycle can be understood via piecewise integration of hyperbolic PDEs, supplemented by shock or simple-wave conditions. Each step—whether at constant pressure or during pressurization/depressurization—admits characteristic equations that are integrable under suitable boundary conditions (e.g., a feed at one end, product at the other, or a closed end).

This equilibrium theory has proven highly useful in analyzing PSA processes because it provides:

• Exact or semi-analytical wave speeds and concentration profiles,
• Simple criteria for identifying shock fronts vs. expansion regions,
• A framework for understanding how separation performance depends on key dimensionless parameters ${\beta }_A$, ${\beta }_B$, and $\beta$.

More detailed calculations can incorporate boundary conditions for each PSA step, track how the shock or wave front moves through the bed, and predict when and where the heavy or light components break through. Although the theory is based on idealized assumptions, it offers valuable insights for process design, performance limits, and optimization of PSA cycles.

Code Overview

The Python implementation follows the equilibrium theory equations described by Knaebel & Hill for simulating wave propagation in a PSA (Pressure Swing Adsorption) column. The core data structures (Region, WaveFront, Characteristic) capture the state of the bed and any traveling fronts or characteristics. The main computation logic resides in the WaveTracker class, which calculates wave speeds at both constant and changing pressures. A higher-level controller, WaveDiagram, orchestrates time-stepping (or pressure-stepping) through multiple PSA steps, updates the column profile, and manages plotting.

Region

A Region stores a constant concentration of the heavy component (y_value) for a segment of the column between z_left and z_right. Multiple Region objects strung together represent the piecewise-constant concentration profile along the column length.

WaveFront

A WaveFront tracks a discontinuity moving through the bed. It has a wave_type that may be "shock" (when a discontinuity steepens) or "simple" (when it spreads out). The attributes z_position, y_ahead, and y_behind indicate the current position of the wavefront and the compositions it separates.

Characteristic

A Characteristic represents a path in the ((z,t)) or ((z,P)) plane corresponding to a particular composition. These paths are stored in arrays (z_values, y_values, t_values, and p_values) for post-processing or plotting. Characteristics are typically generated inside the column and can merge with (or be terminated by) wavefronts.

WaveTracker

The WaveTracker class implements the theoretical equations needed to determine front speeds and characteristic velocities under different operating conditions. It contains methods for:

• Constant-Pressure Wave Speeds. Functions such as get_characteristic_speed_constant_p() and get_shock_speed_constant_p() implement Equations (9)-(10) and (15) of Knaebel & Hill, giving the velocities of simple waves or shocks.
• Changing-Pressure Wave Speeds. Additional methods like get_u_changing_p() and get_shock_speed_changing_p() capture Equations (7) and (16) for situations where the column is undergoing pressurization or depressurization.

By calling these methods, the code can determine how a wavefront or characteristic should move over a small time (or pressure) interval.

WaveDiagram

This class coordinates the overall simulation. It breaks the total PSA process into discrete “steps” (e.g., a constant-pressure feed step, a blowdown step with decreasing pressure, etc.) and handles:

1. Detecting discontinuities in the bed and spawning new wavefronts.
2. Moving wavefronts and characteristics incrementally through time or pressure changes.
3. Rebuilding the piecewise-constant column profile based on where wavefronts have moved.
4. Storing intermediate snapshots and final states for visualization.

When WaveDiagram completes a step, it has an updated column profile plus a record of all wavefronts and characteristics for plotting.

Examples and Usage

Below are simplified examples illustrating how to set up and run PSA steps. The bed length, inlet compositions, pressures, and velocities are all adjustable, making it easy to experiment with different cycle designs.

2-Bed 4-Step Cycle (product for pressurization)

The following code creates a bed of length $L=1.0$ with an initially uniform composition of $y=0.0$. It runs four steps—Pressurization, Feed, Blowdown, and Purge—and then plots the resulting wave diagrams.

Click to expand the code

# Initialize the WaveTracker and WaveDiagram
wave_tracker = WaveTracker(betaA=0.3, betaB=0.7, beta=0.5)
wave_diagram = WaveDiagram(wave_tracker)

# Define the bed and step parameters
L = 1.0
initial_profile = [Region(y_value=0.0, z_left=0.0, z_right=L)]
T_pressurization = 2.0
T_feed = 2.0
T_blowdown = 4.0
T_purge = 2.0
P_low = 1.0
P_high = 5.0
y_feed = 0.8
y_product = 0.0  # pure light

# Pressurization step: increases pressure from P_low to P_high,
# with flow entering from the 'reverse' end (z=L).
wave_diagram.run_changing_pressure_step(
    step_name="Pressurization",
    y_in=y_product,
    p_start=P_low,
    p_end=P_high,
    t_step=T_pressurization,
    z_profile_in=initial_profile,
    flow_dir="reverse",
)

# Feed step: keeps pressure constant, injects heavier component from z=0
wave_diagram.run_constant_pressure_step(
    step_name="Feed",
    u_in=1.8,
    y_in=y_feed,
    t_step=T_feed,
    flow_dir="default",
)

# Blowdown step: decreases pressure from P_high to P_low
wave_diagram.run_changing_pressure_step(
    step_name="Blowdown",
    y_in=None,
    p_start=P_high,
    p_end=P_low,
    t_step=T_blowdown,
    flow_dir="default",
)

# Purge step: at low pressure, flow enters from z=L to desorb heavy component
wave_diagram.run_constant_pressure_step(
    step_name="Purge",
    u_in=2.0,
    y_in=0.0,
    t_step=T_purge,
    flow_dir="reverse",
)

# Generate and display the wave diagram
wave_diagram.plot_steps(["Pressurization","Feed","Blowdown","Purge"])

In this cycle:

1. Pressurization lifts the column from low to high pressure. The code checks whether a wavefront forms at $z=L$ depending on the composition mismatch.
   
2. Feed introduces a heavier component at the inlet $z=0$ under a constant, now high, pressure.
   
3. Blowdown vents the column from the other end, decreasing pressure from $P_{high}$ to $P_{low}$.
   
4. Purge at low pressure pushes a light component backward through the column to remove heavy species.

When the script finishes, the result is a stitched wave diagram showing how composition fronts move step by step.

2-Bed 4-Step Cycle (feed for pressurization)

A similar 4-step sequence can be constructed but this time the feed is used for pressurization. This allows you to see how slight changes in step sequences can lead to very different wavefronts.

4-Bed 9-Step Industrial Cycle

The third example demonstrates a typical industrial multi-step PSA, including multiple depressurizations, equalizations, and repressurizations. The code is longer, but follows the same pattern:

• A series of run_constant_pressure_step and run_changing_pressure_step calls for each sub-step (e.g., Adsorption, partial depressurization, purge).
• Each step specifies relevant inlet compositions, target pressures, and durations.
• After finishing all steps, plot_steps creates a consolidated wave diagram.

In this scenario, you can explore the effects of partial equalizations (PEQ) and depressurizations (DEQ) to see how multi-step cycles shift bed compositions.

How to Run and Adapt

1. Dependencies and Setup. You need Python 3, NumPy, Scipy, and Matplotlib. Save the code as wave_tracking.py (or any name) and run it directly (python wave_tracking.py).
   
2. Modifying the Bed. Adjust or extend initial_profile to reflect different initial loading conditions (multiple Region objects with various y_values).
   
3. Adjusting Equilibrium Parameters. The isotherm parameters kA, kB, and bed porosity epsilon determine the values of betaA, betaB, and beta. Changing these inputs reflects different adsorbent materials.
   
4. Defining Custom Steps. In your own PSA script, call run_constant_pressure_step or run_changing_pressure_step with the times, pressures, flows, and compositions you want. The code automatically spawns shocks or simple waves where discontinuities appear.
   
5. Plotting and Analysis. The plot_steps method creates a stitched diagram of composition fronts across multiple steps. It can be customized or replaced with your own post-processing routines if desired.

By experimenting with different inlet conditions, velocities, and pressure profiles, you can simulate and visualize how composition waves and shocks evolve in a PSA system under various operating strategies.