{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# One-dimensional packed-bed, catalytic-membrane reactor model without streamwise diffusion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The present model simulates heterogeneous catalytic processes inside packed-bed, catalytic membrane reactors. The gas-phase and surface-phase species conservation equations are derived and the system of differential-algebraic equations (DAE) is solved using the scikits.odes.dae IDA solver."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Runnning Cantera version: 2.6.0\n"
]
}
],
"source": [
"# Import Cantera and scikits\n",
"import numpy as np\n",
"from scikits.odes import dae\n",
"import cantera as ct\n",
"import matplotlib.pyplot as plt\n",
"\n",
"%matplotlib inline\n",
"%config InlineBackend.figure_formats = [\"svg\"]\n",
"print(\"Runnning Cantera version: \" + ct.__version__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Methodology\n",
"\n",
"One-dimensional, steady-state catalytic-membrane reactor model with surface chemistry is developed to analyze species profiles along the length of a packed-bed, catalytic membrane reactor. The same model can further be simplified to simulate a simple packed-bed reactor by excluding the membrane. The example here demonstrates the one-dimensional reactor model explained by G. Kogekar, et al., \"Efficient and robust computational models of heterogeneous catalysis in catalytic membrane reactors,\" (2022), In preparation [1]."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Governing equations\n",
"\n",
"Assuming steady-state, one-dimensional flow within the packed bed, total-mass, species mass and energy conservation may be stated as [2]\n",
"\n",
"$$\n",
" \\frac{d(\\rho u)}{dz} = \\sum_{k=1}^{K_{\\mathrm{g}}} \\dot {s}_k W_k A_{\\mathrm{s}} + \\frac{P_{\\mathrm b}}{A_{\\mathrm b}} j_{k_{\\mathrm M}}, \\\\\n",
" \\rho u \\frac{dY_k}{dz} + A_{\\mathrm{s}} Y_k \\sum_{k=1}^{K_{\\mathrm{g}}} \\dot {s}_k W_k = A_{\\mathrm s} \\dot {s}_k W_k + \\delta_{k, k_{\\mathrm M}} \\frac{P_{\\mathrm b}}{A_{\\mathrm b}} j_{k_{\\mathrm M}}, \\\\\n",
" \\rho u c_{\\mathrm p} \\frac{dT}{dz} + \\sum_{k=1}^{K_{\\mathrm g}} h_k (\\phi_{\\mathrm g} \\dot {\\omega}_k + A_{\\mathrm s} \\dot {s}_k) W_k = \\hat h \\frac{P_{\\mathrm b}}{A_{\\mathrm b}}(T_{\\mathrm w} - T) + \\delta_{k, k_{\\mathrm M}} \\frac{P_{\\mathrm b}}{A_{\\mathrm b}} h_{k_{\\rm M}} j_{k_{\\mathrm M}}.\n",
"$$\n",
"\n",
"The fractional coverages of the $K_{\\rm s}$ surface adsorbates $\\theta_k$ must satisfy \n",
"$$\n",
" \\dot {s}_k = 0, {\\ \\ \\ \\ \\ \\ } (k = 1,\\ldots, K_{\\mathrm s}),\n",
"$$\n",
"which, at steady state, requires no net production/consumption of surface species by the heterogeneous reactions [3].\n",
"\n",
"The pressure within the bed is calculated as\n",
"$$ \n",
" \\frac{dp}{dz} = - \\left( \\frac{\\phi_{\\mathrm{g}} \\mu}{\\beta_{\\mathrm{g}}} \\right) u,\n",
"$$\n",
"where $\\mu$ is the gas-phase mixture viscosity. The packed-bed permeability $\\beta_{\\rm g}$ is evaluated using the Kozeny-Carman relationship as \n",
"$$\n",
" \\beta_{\\mathrm g} = \\frac{\\phi_{\\mathrm{g}}^3 D_{\\mathrm p}^2}{72 \\tau_{\\mathrm{g}}(1 - \\phi_{\\mathrm{g}})^2},\n",
"$$\n",
"where $\\phi_{\\mathrm{g}}$, $\\tau_{\\mathrm{g}}$, and $D_{\\mathrm {p}}$ are the bed porosity, tortuosity, and particle diameter, respectively.\n",
"\n",
"The independent variable in these conservation equations is the position $z$ along the reactor length. The dependent variables include total mass flux $\\rho u$, pressure $p$, temperature $T$, gas-phase mass fractions $Y_k$, and surfaces coverages $\\theta_k$. Gas-phase fluxes through the membrane are represented as $j_{k, {\\mathrm M}}$. Geometric parameters $A_{\\mathrm s}$, $P_{\\mathrm b}$, and $A_{\\mathrm b}$ represent the catalyst specific surface area (dimensions of surface area per unit volume), reactor perimeter, and reactor cross-sectional flow area, respectively. Other parameters include bed porosity $\\phi_{\\mathrm g}$ and gas-phase species molecular weights $W_k$. The gas density $\\rho$ is evaluated using the equation of state (ideal Eos, RK or PR EoS).\n",
"\n",
"If a perm-selective membrane is present, then $j_{k_{\\mathrm M}}$ represents the gas-phase flux through the membrane and ${k_{\\mathrm M}}$ is the gas-phase species that permeates through the membrane. The Kronecker delta, $\\delta_{k, k_{\\mathrm M}} = 1$ for the membrane-permeable species and $\\delta_{k, k_{\\mathrm M}} = 0$ otherwise. The membrane flux is calculated using Sievert's law as\n",
"$$\n",
" {j_{k_{\\mathrm M}}}^{\\text{Mem}} = \\frac{B_{k_{\\mathrm M}}}{t} \\left ( p_{k_{\\mathrm M} {\\text{, mem}}}^\\alpha - p_{k_{\\mathrm M} \\text{, sweep}}^\\alpha \\right ) W_{k_{\\rm M}}\n",
"$$\n",
"where $B_{k_{\\mathrm M}}$ is the membrane permeability, $t$ is the membrane thickness. $p_{k_{\\mathrm M} \\text{, mem}}$ and $p_{k_{\\mathrm M} \\text{, sweep}}$ represent perm-selective species partial pressures within the packed-bed and the exterior sweep channel. The present model takes the pressure exponent $\\alpha$ to be unity. The membrane flux for all other species ($ k \\neq k_{\\mathrm M}$) is zero.\n",
"\n",
"### Chemistry mechanism\n",
"This example uses a detailed 12-step elementary micro-kinetic reaction mechanism that describes ammonia formation and decomposition kinetics over the Ru/Ba-YSZ catalyst. The reaction mechanism is developed and validated using measured performance in a laboratory-scale packed-bed reactor [4]. This example also incorporates a Pd-based H2 perm-selective membrane.\n",
"\n",
"\n",
"### Solver\n",
"\n",
"Above governing equations represent a complete solution for a steady-state packed-bed, catalytic membrane reactor model. The dependent variables are the mass-flux $\\rho u$, species mass-fractions $Y_k$, pressure $p$, temperature $T$, and surface coverages $\\theta_k$. The equation of state is used to obtain the mass density, $\\rho$.\n",
"\n",
"The governing equations form an initial value problem (IVP) in a differential-algebraic (DAE) form as follows:\n",
"$$\n",
" f(z,{\\bf{y}}, \\bf {y'}, c) = 0,\n",
"$$\n",
"where $\\bf{y}$ and $\\bf{y'}$ represent the solution vector and its derivative vector, respectively. All other constants such as reference temperature, chemical constants, etc. are incorporated in vector $c$ (Refer to [5] for more details). This type of DAE system in this example is solved using the `scikits.odes.dae` IDA solver.\n",
"> 1. G. Kogekar, H. Zhu, R.J. Kee, 'Efficient and robust computational models of heterogeneous catalysis in catalytic membrane reactors', (2022), In preparation\n",
"> 2. B. Kee, C. Karakaya, H. Zhu, S. DeCaluwe, and R.J. Kee, 'The Influence of Hydrogen-Permeable Membranes and Pressure on Methane Dehydroaromatization in Packed-Bed Catalytic Reactors', Industrial & Engineering Chemistry Research (2017) 56, 13:3551 - 3559\n",
"> 3. R.J. Kee, M.E. Coltrin, P. Glarborg, and H. Zhu, 'Chemically Reacting Flow: Theory, Modeling and Simulation', Wiley (2018)\n",
"> 4. Z. Zhang, C. Karakaya, R.J. Kee, J. Douglas Way, C. Wolden,, 'Barium-Promoted Ruthenium Catalysts on Yittria-Stabilized Zirconia Supports for Ammonia Synthesis', ACS Sustainable Chemistry & Engineering (2019) 7:18038 - 18047\n",
"> 5. G. Kogekar, 'Computationally efficient and robust models of non-ideal thermodynamics, gas-phase kinetics and heterogeneous catalysis in chemical reactors' (2021), Doctoral dissertation. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define gas-phase and surface-phase species"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Number of gas-phase species = 4\n",
"Number of surface-phase species = 6\n",
"Number of variables = 13\n"
]
}
],
"source": [
"# Import the reaction mechanism for Ammonia synthesis/decomposition on Ru-Ba/YSZ catalyst\n",
"mechfile = \"data/Ammonia-Ru-Ba-YSZ.yaml\"\n",
"# Import the models for gas-phase\n",
"gas = ct.Solution(mechfile, \"gas\")\n",
"# Import the model for surface-phase\n",
"surf = ct.Interface(mechfile, \"Ru_surface\", [gas])\n",
"\n",
"# Other parameters\n",
"n_gas = gas.n_species # number of gas species\n",
"n_surf = surf.n_species # number of surface species\n",
"n_gas_reactions = gas.n_reactions # number of gas-phase reactions\n",
"\n",
"# Set offsets of dependent variables in the solution vector\n",
"offset_rhou = 0\n",
"offset_p = 1\n",
"offset_T = 2\n",
"offset_Y = 3\n",
"offset_Z = offset_Y + n_gas\n",
"n_var = offset_Z + n_surf # total number of variables (rhou, P, T, Yk and Zk)\n",
"\n",
"print(\"Number of gas-phase species = \", n_gas)\n",
"print(\"Number of surface-phase species = \", n_surf)\n",
"print(\"Number of variables = \", n_var)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define reactor geometry and operating conditions"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Modeling packed-bed, catalytic-membrane reactor...\n",
"H2 permeable membrane is present.\n"
]
}
],
"source": [
"# Reactor geometry\n",
"L = 5e-2 # length of the reactor (m)\n",
"R = 5e-3 # radius of the reactor channel (m)\n",
"phi = 0.5 # porosity of the bed (-)\n",
"tau = 2.0 # tortuosity of the bed (-)\n",
"dp = 3.37e-4 # particle diameter (m)\n",
"As = 3.5e6 # specific surface area (1/m)\n",
"\n",
"# Energy (adiabatic or isothermal)\n",
"solve_energy = True # True: Adiabatic, False: isothermal\n",
"\n",
"# Membrane (True: membrane, False: no membrane)\n",
"membrane_present = True\n",
"membrane_perm = 1e-15 # membrane permeability (kmol*m3/s/Pa)\n",
"thickness = 3e-6 # membrane thickness (m)\n",
"membrane_sp_name = \"H2\" # membrane-permeable species name\n",
"p_sweep = 1e5 # partial pressure of permeable species in the sweep channel (Pa)\n",
"permeance = membrane_perm / thickness # permeance of the membrane (kmol*m2/s/Pa)\n",
"\n",
"if membrane_present:\n",
" print(\"Modeling packed-bed, catalytic-membrane reactor...\")\n",
" print(membrane_sp_name, \"permeable membrane is present.\")\n",
"\n",
"# Get required properties based on the geometry and mechanism\n",
"W_g = gas.molecular_weights # vector of molecular weight of gas species\n",
"vol_ratio = phi / (1 - phi)\n",
"eff_factor = phi / tau # effective factor for permeability calculation\n",
"# permeability based on Kozeny-Carman equation\n",
"B_g = B_g = vol_ratio**2 * dp**2 * eff_factor / 72 \n",
"area2vol = 2 / R # area to volume ratio assuming a cylindrical reactor\n",
"D_h = 2 * R # hydraulic diameter\n",
"membrane_sp_ind = gas.species_index(membrane_sp_name)\n",
"\n",
"# Inlet operating conditions\n",
"T_in = 673 # inlet temperature [K]\n",
"p_in = 5e5 # inlet pressure [Pa]\n",
"v_in = 0.001 # inlet velocity [m/s]\n",
"T_wall = 723 # wall temperature [K]\n",
"h_coeff = 1e2 # heat transfer coefficient [W/m2/K]\n",
"\n",
"# Set gas and surface states\n",
"gas.TPX = T_in, p_in, \"NH3:0.99, AR:0.01\" # inlet composition\n",
"surf.TP = T_in, p_in\n",
"Yk_0 = gas.Y\n",
"rhou0 = gas.density * v_in\n",
"\n",
"# Initial surface coverages\n",
"# advancing coverages over a long period of time to get the steady state.\n",
"surf.advance_coverages(1e10) \n",
"Zk_0 = surf.coverages"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Define residual function required for IDA solver"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def residual(z, y, yPrime, res):\n",
" \"\"\"Solution vector for the model\n",
" y = [rho*u, p, T, Yk, Zk]\n",
" yPrime = [d(rho*u)dz, dpdz, dTdz, dYkdz, dZkdz]\n",
" \"\"\"\n",
" # Get current thermodynamic state from solution vector and save it to local variables.\n",
" rhou = y[offset_rhou] # mass flux (density * velocity)\n",
" Y = y[offset_Y : offset_Y + n_gas] # vector of mass fractions\n",
" Z = y[offset_Z : offset_Z + n_surf] # vector of site fractions \n",
" p = y[offset_p] # pressure\n",
" T = y[offset_T] # temperature\n",
"\n",
" # Get derivatives of dependent variables\n",
" drhoudz = yPrime[offset_rhou] \n",
" dYdz = yPrime[offset_Y : offset_Y + n_gas] \n",
" dZdz = yPrime[offset_Z : offset_Z + n_surf] \n",
" dpdz = yPrime[offset_p] \n",
" dTdz = yPrime[offset_T] \n",
"\n",
" # Set current thermodynamic state for the gas and surface phases\n",
" # Note: use unnormalized mass fractions and site fractions to avoid over-constraining the system\n",
" gas.set_unnormalized_mass_fractions(Y)\n",
" gas.TP = T, p\n",
" surf.set_unnormalized_coverages(Z)\n",
" surf.TP = T, p\n",
"\n",
" # Calculate required variables based on the current state\n",
" coverages = surf.coverages # surface site coverages\n",
" # heterogeneous production rate of gas species\n",
" sdot_g = surf.get_net_production_rates(\"gas\")\n",
" # heterogeneous production rate of surface species\n",
" sdot_s = surf.get_net_production_rates(\"Ru_surface\")\n",
" wdot_g = np.zeros(n_gas)\n",
" # specific heat of the mixture\n",
" cp = gas.cp_mass \n",
" # partial enthalpies of gas species\n",
" hk_g = gas.partial_molar_enthalpies \n",
"\n",
" if n_gas_reactions > 0:\n",
" # homogeneous production rate of gas species\n",
" wdot_g = gas.net_production_rates \n",
" mu = gas.viscosity # viscosity of the gas-phase\n",
"\n",
" # Calculate density using equation of state\n",
" rho = gas.density\n",
"\n",
" # Calculate flux term through the membrane\n",
" # partial pressure of membrane-permeable species\n",
" memsp_pres = (p * gas.X[membrane_sp_ind]) \n",
" # negative sign indicates the flux going out\n",
" membrane_flux = (-permeance * (memsp_pres - p_sweep) / W_g[membrane_sp_ind]) \n",
"\n",
" # Conservation of total-mass\n",
" # temporary variable\n",
" sum_continuity = As * np.sum(sdot_g * W_g) + phi * np.sum(wdot_g * W_g) \n",
" res[offset_rhou] = (drhoudz - sum_continuity - area2vol * membrane_flux * membrane_present)\n",
"\n",
" # Conservation of gas-phase species \n",
" res[offset_Y:offset_Y+ n_gas] = (dYdz + (Y * sum_continuity - phi * np.multiply(wdot_g,W_g) \n",
" - As * np.multiply(sdot_g,W_g)) / rhou)\n",
" res[offset_Y + membrane_sp_ind] -= area2vol * membrane_flux * membrane_present\n",
"\n",
" # Conservation of site fractions (algebraic contraints in this example)\n",
" res[offset_Z : offset_Z + n_surf] = sdot_s\n",
"\n",
" # For the species with largest site coverage (k_large), solve the constraint equation: sum(Zk) = 1\n",
" # The residual function for 'k_large' would be 'res[k_large] = 1 - sum(Zk)'\n",
" # Note that here sum(Zk) will be the sum of coverages for all surface species, including the 'k_large' species.\n",
" ind_large = np.argmax(coverages)\n",
" res[offset_Z + ind_large] = 1 - np.sum(coverages)\n",
"\n",
" # Conservation of momentum\n",
" u = rhou / rho\n",
" res[offset_p] = dpdz + phi * mu * u / B_g\n",
"\n",
" # Conservation of energy\n",
" res[offset_T] = dTdz - 0 # isothermal condition\n",
" # Note: One can just not solve the energy equation by keeping temperature constant. \n",
" # But since 'T' is used as the dependent variable, the residual is res[T] = dTdz - 0 \n",
" # One can also write res[T] = 0 directly, but that leads to a solver failure due to singular jacobian\n",
"\n",
" if solve_energy:\n",
" conv_term = (4 / D_h) * h_coeff * (T_wall - T) * (2 * np.pi * R)\n",
" chem_term = np.sum(hk_g * (phi * wdot_g + As * sdot_g))\n",
" res[offset_T] -= (conv_term - chem_term) / (rhou * cp)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Calculate the spatial derivatives at the inlet that will be used as the initial conditions for the IDA solver\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Initialize yPrime to 0 and call residual to get initial derivatives\n",
"y0 = np.hstack((rhou0, p_in, T_in, Yk_0, Zk_0))\n",
"yprime0 = np.zeros(n_var)\n",
"res = np.zeros(n_var)\n",
"residual(0, y0, yprime0, res)\n",
"yprime0 = -res"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Solve the system of DAEs using ida solver"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SolverReturn(flag=, values=SolverVariables(t=0.05008188075110137, y=array([1.40862546e-03, 4.99997977e+05, 6.85053268e+02, 8.87695903e-02,\n",
" 4.76156874e-01, 4.11794146e-01, 2.31457254e-02, 3.31389994e-03,\n",
" 9.93036016e-01, 2.72895231e-03, 8.58521045e-04, 3.67135998e-08,\n",
" 6.25735914e-05]), ydot=array([-1.00712547e-02, -4.63384892e+01, 2.35197310e+01, 1.93773529e+00,\n",
" -1.09694323e+01, 9.02163101e+00, -4.93684753e-06, 1.65140796e-01,\n",
" -3.26729842e-01, 1.54845328e-01, 5.55655962e-03, 4.81233014e-07,\n",
" 1.18667692e-03])), errors=SolverVariables(t=None, y=None, ydot=None), roots=SolverVariables(t=None, y=None, ydot=None), tstop=SolverVariables(t=None, y=None, ydot=None), message='Successful function return.')\n"
]
}
],
"source": [
"solver = dae(\n",
" \"ida\",\n",
" residual,\n",
" first_step_size=1e-15,\n",
" atol=1e-14, # absolute tolerance for solution\n",
" rtol=1e-06, # relative tolerance for solution\n",
" algebraic_vars_idx=[np.arange(offset_Y + n_gas, offset_Z + n_surf, 1)],\n",
" max_steps=8000,\n",
" one_step_compute=True,\n",
" old_api=False, # forces use of new api (namedtuple)\n",
")\n",
"\n",
"distance = []\n",
"solution = []\n",
"state = solver.init_step(0.0, y0, yprime0)\n",
"\n",
"# Note that here the variable t is an internal varible used in scikits. In this example, it represents\n",
"# the natural variable z, which corresponds to the axial distance inside the reactor.\n",
"while state.values.t < L:\n",
" distance.append(state.values.t)\n",
" solution.append(state.values.y)\n",
" state = solver.step(L)\n",
"\n",
"distance = np.array(distance)\n",
"solution = np.array(solution)\n",
"print(state)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Plot results"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
"