Isothermal design of batch, CSTR, PFR and packed-bed reactors

With a rate law and a stoichiometric table, the design equations of lessons 1 and 2 can be solved for any isothermal reactor. Fogler's textbook calls this the algorithm of reaction engineering: write the mole balance, the rate law and the stoichiometry, combine them, and evaluate. This lesson applies it to two plant problems. The first is a batch: Dev Raman's laboratory wants to know how long its one-liter esterification runs will take, and the operations group asks whether a batch reactor could make ester while the continuous line is being commissioned. The second is the ethylene pilot's packed bed, where pressure drop through the catalyst competes with the reaction.

The algorithm

When the combined equation has an analytic integral, it gives a formula; when it does not, as soon as pressure drop or several reactions appear, it is integrated numerically. The code in this course uses the classical fourth-order Runge-Kutta method, which takes four slope evaluations per step and is accurate enough for reactor design with modest step sizes.

Batch esterification

For the equimolar esterification at constant volume, the batch balance is \(dX/dt = kC_{A0}[(1 - X)^2 - X^2/K]\), with \(kC_{A0} = 1.3\times10^{-3}\) per second for a well-stirred mixture with the same amount of resin per liter as the pre-reactor bed. The equation has an analytic solution, but the numerical one is easier to adapt, so the code integrates it.

The reaction reaches 30 % in under 6 minutes and 50 % in 14 minutes, then slows as it approaches equilibrium: 60 % takes 25 minutes and 65 % takes 42. A batch reactor in production would also need time to charge the reactants, heat them to 80 °C, filter off the resin or let it settle, cool and empty the vessel, perhaps an hour or more. Dev's laboratory runs will take about half an hour of reaction time each, and the operations group's interim batch scheme would make one batch every two hours or so.

Packed bed with pressure drop

The ethylene pilot unit passes ethanol vapor diluted with steam, 40 % ethanol at 3 bar inlet pressure and 350 °C, through a bed of solid acid catalyst. At these conditions the dehydration is effectively irreversible and first order in ethanol, so per kilogram of catalyst \(-r_A' = k'C_A\), with \(k' = 1.8\times10^{-4}\) m³/(kg s), an illustrative value for the pilot catalyst. Ethanol feed is 4.25 kmol/h, about 1,000 metric tons of ethylene a year if all converted.

Gas flowing through a packed bed loses pressure by friction. The Ergun equation describes this loss in terms of particle size, void fraction, gas density and velocity. For an isothermal bed it reduces to a simple form for the pressure ratio \(p = P/P_0\):

where \(\alpha\) collects the bed properties. The \((1 + \varepsilon X)\) term appears because the expanding gas flows faster. Pressure drop matters to the reaction because the concentration is proportional to pressure: \(C_A = C_{A0}(1 - X)p/(1 + \varepsilon X)\). The mole balance and the Ergun equation must be integrated together.

# Isothermal design: batch esterification time, and a packed bed for ethanol dehydration with Ergun pressure drop
import math
# 1. Batch esterification, 1:1, K = 4, k CA0 = 1.3e-3 1/s at 80 C (illustrative, from lesson 2)
kc = 1.5e-4 * 8.7
def dXdt(X): return kc * ((1 - X) ** 2 - X ** 2 / 4.0)
t, X, h = 0.0, 0.0, 1.0
marks = [0.3, 0.5, 0.6, 0.65]
while marks:
    k1 = dXdt(X); k2 = dXdt(X + h * k1 / 2); k3 = dXdt(X + h * k2 / 2); k4 = dXdt(X + h * k3)
    X += h * (k1 + 2 * k2 + 2 * k3 + k4) / 6; t += h
    if X >= marks[0]:
        print(f"batch reaches X = {marks.pop(0):.2f} after {t / 60:5.1f} min"); 
print("a full batch cycle adds charging, heating, cooling and emptying time")

# 2. Packed bed: C2H5OH -> C2H4 + H2O, first order in ethanol, isothermal at 350 C
R = 0.08314
P0, T = 3.0, 623.15          # bar inlet, K
yA0 = 0.4; eps = yA0           # 40 % ethanol in steam
FT0 = 10.6 / 3600              # kmol/s total feed (4.25 kmol/h ethanol)
FA0 = yA0 * FT0
CA0 = yA0 * P0 / (R * T)       # kmol/m3
kp = 1.8e-4                    # m3/(kg cat s), illustrative first-order constant at 350 C
alpha = 2.0e-4                 # 1/kg, Ergun pressure-drop parameter for the chosen particles (illustrative)

def f(W, X, y, dp=True):
    CA = CA0 * (1 - X) / (1 + eps * X) * y
    dX = kp * CA / FA0
    dy = -alpha / (2 * y) * (1 + eps * X) if dp else 0.0
    return dX, dy

def run(Wmax, dp):
    W, X, y, h = 0.0, 0.0, 1.0, 1.0
    out = {}
    while W < Wmax - 1e-9:
        a = f(W, X, y, dp); b = f(W + h / 2, X + h * a[0] / 2, y + h * a[1] / 2, dp)
        c = f(W + h / 2, X + h * b[0] / 2, y + h * b[1] / 2, dp); d = f(W + h, X + h * c[0], y + h * c[1], dp)
        X += h * (a[0] + 2 * b[0] + 2 * c[0] + d[0]) / 6; y += h * (a[1] + 2 * b[1] + 2 * c[1] + d[1]) / 6; W += h
        if y < 0.05: break
        if int(W) % 500 == 0: out[int(W)] = (X, y)
    return out
noDP = run(2500, False); DP = run(2500, True)
print("W (kg)   X without dP   X with dP   P/P0")
for W in (500, 1000, 1500, 2000, 2500):
    print(f"{W:6d}      {noDP[W][0]:.3f}        {DP[W][0]:.3f}      {DP[W][1]:.2f}")
# catalyst needed for X = 0.98 without pressure drop (analytic): W = FA0/(k CA0) [(1+eps) ln(1/(1-X)) - eps X]
X = 0.98
Wa = FA0 / (kp * CA0) * ((1 + eps) * math.log(1 / (1 - X)) - eps * X)
print(f"analytic catalyst mass for X = 0.98 with no pressure drop: {Wa:.0f} kg")

batch reaches X = 0.30 after   5.6 min
batch reaches X = 0.50 after  14.0 min
batch reaches X = 0.60 after  24.9 min
batch reaches X = 0.65 after  42.1 min
a full batch cycle adds charging, heating, cooling and emptying time
W (kg)   X without dP   X with dP   P/P0
   500      0.774        0.764      0.94
  1000      0.939        0.928      0.86
  1500      0.983        0.975      0.78
  2000      0.995        0.990      0.68
  2500      0.999        0.996      0.57
analytic catalyst mass for X = 0.98 with no pressure drop: 1436 kg

Reading the bed results

Without pressure drop, 1,000 kg of catalyst gives 93.9 % conversion and 1,436 kg reaches 98 %, matching the analytic formula \(W = (F_{A0}/k'C_{A0})[(1 + \varepsilon)\ln(1/(1-X)) - \varepsilon X]\). With pressure drop, conversion lags slightly, 92.8 % at 1,000 kg, and the pressure falls to 86 % of the inlet value; at 2,000 kg it is down to 68 %. The conversion penalty is small because the reaction is first order and fast. The larger cost is downstream: ethylene leaves at lower pressure and must be compressed further before separation and storage.

Choosing particles is a trade-off. Smaller particles give more external surface and less internal diffusion resistance, which raises the effective rate, but they raise the pressure drop roughly in inverse proportion to particle size, and more steeply at high flow. Engineers often use shaped extrudates, rings or trilobes, which give a short diffusion path with an open bed. A wider, shorter bed lowers the velocity and the pressure drop for the same catalyst mass, at the cost of a larger vessel and a harder job distributing the gas evenly across it.

Isothermal CSTR and PFR formulas

For reference, the simplest cases have closed forms. For a first-order liquid reaction with \(k\tau = Da\), the Damköhler number, a CSTR gives \(X = Da/(1 + Da)\) and a PFR gives \(X = 1 - e^{-Da}\). For a second-order liquid reaction with equal feeds and \(Da = kC_{A0}\tau\), a CSTR gives \(X = [(1 + 2Da) - \sqrt{1 + 4Da}]/(2Da)\) and a PFR gives \(X = Da/(1 + Da)\). The Damköhler number compares the reaction rate with the flow rate through the reactor; values much less than 0.1 mean little conversion, values above 10 mean nearly complete conversion of an irreversible reaction.

Two first-order CSTRs in series of equal size each with \(Da\) give \(X = 1 - 1/(1 + Da)^2\); \(n\) tanks give \(1 - 1/(1 + Da)^n\). As \(n\) grows with the total space time fixed, this approaches \(1 - e^{-Da_{\text{total}}}\), the PFR result, which is the analytic version of the Levenspiel staircase in lesson 2.

Scale-up of a batch reaction

Scaling a batch reaction from Dev's one-liter flask to a production vessel keeps the kinetics but changes almost everything else. Heat transfer area per unit volume falls as the vessel grows, roughly in proportion to one over the diameter, so a reaction that held its temperature easily in a flask may need a cooling coil or an external exchanger in a large tank. Mixing times grow too. For the nearly thermoneutral esterification these effects are minor; for the strongly exothermic reactions of lesson 9 they decide whether a process is safe. Every scale-up therefore starts by checking that the rate measured in the laboratory was not limited by mixing or heat transfer, and that the production vessel can remove the heat at the same rate the reaction will release it.

Checking a numerical answer

A numerical integration is only as good as its step size and its equations. Theo checks the code three ways before he reports anything. First, he turns off the pressure drop and compares the result with the analytic formula: the numerical 98.3 % at 1,500 kg and the formula's 1,436 kg for 98 % agree. Second, he halves the step size and confirms that the answers do not change in the third decimal place. Third, he checks the limits: with no catalyst there is no conversion, with very large catalyst mass the conversion approaches 1, and the pressure ratio never rises along the bed. Ingrid asks for this kind of check on every calculation that goes into a design basis, because a single wrong sign in a differential equation can give a plausible-looking but wrong number, and nobody downstream will rerun the code.

The same habits apply to simulator results. A process simulator's reactor model is a set of differential equations solved numerically with the user's kinetics. Checking a simulated reactor against a simple hand calculation at one condition is the fastest way to catch a unit error in a rate constant or a wrong basis for the rate, the two most common mistakes in practice.