PFR Stationary 1 dimension

Plug flow reactor in steady state operation.

import numpy as np

import matplotlib.pyplot as plt

import reactord as rd

from scipy.constants import R

Operating the PFR

Example: Adiabatic Liquid-Phase Isomerization of Normal Butane

Fogler 6th ed. Example 11.3.

Reaction:

\[Butane + iPentane \leftrightarrow iButane + iPentane\]

Kinetic law:

\[r = k_dC_{But} - k_iC_{i-But}\]

We want to process 163 kmol/h and achieve 70% conversion of n-butane from a mixture 90 mol % n-butane and 10 mol % i-pentane, which is considered an inert. The feed enters at 330 K

Specific reaction rate at 360K = 31.1 \(\ {h^{-1}}\)

Kinetic constant at 60 °C = 3.03

Enthalpy reaction (dH) = -6900 \(\frac {J} {mol n-butane}\)

Activation energy (e)= 65.7 \(\frac {kJ} {mol}\)

\(C_{p-nBut} = C_{p-iBut}\) = 141 \(\frac {J} {mol·K}\)

\(C_{p-iPen}\) = 161 \(\frac {J} {mol·K}\)

\(C_{A0}=9.3\) \(\frac {kmol} {mol^3}\)

Reactor: liquid phase, Adiabatic and isobaric PFR Pressure not specified. It is assumed as 101325 Pa (pressure will not be used)

For convenience the reactor will be import from the modules as pfr

import reactord.flowreactors.stationary_1d.pfr as pfr

Substance definitions

For a liquid phase, adiabatic and isobaric PFR the volume and cp of the pure substances are necessary to evaluate the balances.

def volume(temperature, pressure):
    f_mol = 163 * 1000 / 3600  # mol / s
    f_volumetric = 100_000 * 0.00378541 / 24 / 3600  # m3 / s

    rho = f_mol / f_volumetric
    v = 1 / rho  # m3 / mol
    return np.full(np.size(temperature), v)


def cp_butane(temperature, pressure):
    return np.full(np.size(temperature), 141)  # J / mol / K


def cp_pentane(temperature, pressure):
    return np.full(np.size(temperature), 161)  # J / mol / K

but = rd.Substance("but", volume_liquid=volume, heat_capacity_liquid=cp_butane)

ibut = rd.Substance(
    "i-but", volume_liquid=volume, heat_capacity_liquid=cp_butane
)

ipen = rd.Substance(
    "i-pen", volume_liquid=volume, heat_capacity_liquid=cp_pentane
)

Mixture definition

mix = rd.mix.IdealSolution([but, ibut, ipen])

Kinetic definition

# Reaction enthalpy
dh = -6900  # J / mol


# Reaction rate
def r_rate(c, t, cons):
    k360, e, keq60, dh = cons["k360"], cons["e"], cons["keq60"], cons["dh"]

    kd_t = k360 * np.exp(e / R * (1 / 360 - 1 / t))
    keq_t = keq60 * np.exp(dh / R * (1 / (60 + 273.15) - 1 / t))

    rd = kd_t * c["but"]
    ri = kd_t / keq_t * c["i-but"]

    return rd - ri


# Kinetic object
kinetic = rd.Kinetic(
    mix=mix,
    reactions={"r1": {"eq": but > ibut, "rate": r_rate, "DH": dh}},
    kinetic_constants={
        "k360": 31.1 / 3600,
        "e": 65.7 * 1000,
        "keq60": 3.03,
        "dh": dh,
    },
    rates_argument="concentration",
)

For each reactor, different mass, energy and pressure balances can be chosen.

A reactor may be operated in terms of its catalyst operation, thermal operation, and pressure operation.

Catalyst operation	Thermal opration	Pressure operation
Homogeneous	Isothermal	Isobaric
Heterogeneous	Adiabatic	Non-isobaric
	Non-isothermal

For the piston flow reactor, the list of implemented operations is:

• Homogeneous

  • Isothermal

    • Isobaric

    • Non-isobaric

  • Adiabatic

    • Isobaric

    • Non-isobaric

Mass Balance Set

f_mol = 163 * 1000 / 3600  # mol / s

mb = pfr.mass_balances.MolarFlow(
    molar_flows_in={"but": f_mol * 0.9, "i-but": 0, "i-pen": f_mol * 0.1},
)

Energy Balance Set

eb = pfr.energy_balances.Adiabatic(temperature_in_or_out={"in": 330})

Pressure Balance Set

pb = pfr.pressure_balances.Isobaric(101325)

CHECK THE OTHER POSIBBLE MASS, ENERGY AND PRESSURE BALANCES IN THE ReactorD API’s documentation

Set and Simulate the reactor

reactor = pfr.PFR(
    kinetic=kinetic,
    reactor_length=5,
    transversal_area=1,
    grid_size=100,
    mass_balance=mb,
    energy_balance=eb,
    pressure_balance=pb,
)

reactor.simulate(1e-5)

Simulations results are stored in a pandas data frame.

reactor.sim_df

	z	but	i-but	i-pen	temperature	pressure
0	0.000000	40.750000	0.000000	4.527778	330.000000	101325.0
1	0.050505	40.191385	0.558615	4.527778	330.595306	101325.0
2	0.101010	39.618720	1.131280	4.527778	331.205585	101325.0
3	0.151515	39.031625	1.718375	4.527778	331.831243	101325.0
4	0.202020	38.429737	2.320263	4.527778	332.472664	101325.0
...	...	...	...	...	...	...
226	4.932660	11.643340	29.106660	4.527778	361.018466	101325.0
227	4.949495	11.643324	29.106676	4.527778	361.018483	101325.0
228	4.966330	11.643310	29.106690	4.527778	361.018499	101325.0
229	4.983165	11.643295	29.106705	4.527778	361.018514	101325.0
230	5.000000	11.643282	29.106718	4.527778	361.018528	101325.0

231 rows × 6 columns

z = reactor longitude [m]

but = molar flow of butane [mol/s]

i-but = molar flow of iso-butane [mol/s]

i-pen = molar flow of iso-pentane [mol/s]

temperature = reactor’s temperature [K]

pressure = reactor’s pressure [Pa]

Ploting results

results = reactor.sim_df
rd_fbut = results["but"]
rd_x = (f_mol * 0.9 - rd_fbut) / (f_mol * 0.9)

fig1, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(20, 5))

plt.rcParams["font.size"] = 15

ax1.set_xlabel("Reactor Volume [$m^{3}$]")
ax1.set_ylabel("Concentrations (mol/$m^{3}$)")
ax1.set_title("Concentrations")
ax1.plot(results.z, results["but"], "-k", label="butane", linewidth=2.5)
ax1.plot(results.z, results["i-but"], "-g", label="iso-butane", linewidth=2.5)

ax1.plot(
    results.z, results["i-pen"], "--c", label="iso-pentane", linewidth=1.8
)

ax1.legend()

ax2.set_xlabel("Reactor Volume [$m^{3}$]")
ax2.set_ylabel("Temperature (K)")
ax2.set_title("Temperature")
ax2.plot(results.z, results.temperature, "-r")

ax3.set_xlabel("Reactor Volume [$m^{3}$]")
ax3.set_ylabel("Pressure (K)")
ax3.set_title("Pressure")
ax3.plot(
    results.z,
    results.pressure,
    "-b",
);

The reactor object can print the LaTex representation for the equation that solves and in the dunder method repr the LaTex code to be used.

reactor.irepr
reactor

Mass balance:

$\displaystyle \frac{1}{a_t}\frac{dF_i}{dz}=r_i$

Reactor and refrigerant energy balances:

$\displaystyle \frac{dT}{dz}=\frac{\sum\Delta{H_j}_{(TP)}r_{j}}{{c_p}_{mix}{\sum}F_i}$

$\displaystyle \frac{dT_r}{dz}=0$

Pressure balance:

$\displaystyle \frac{dP}{dz}=0$

\frac{1}{a_t}\frac{dF_i}{dz}=r_i
\frac{dT}{dz}=\frac{\sum\Delta{H_j}_{(TP)}r_{j}}{{c_p}_{mix}{\sum}F_i}
\frac{dT_r}{dz}=0
\frac{dP}{dz}=0