Cantera/Reactor Equations
From charlesreid1
The thermochemical state of a gas inside a gas reactor is set by specifying the gas's energy state, pressure, and composition. Cantera does this by solving equations for a the gas internal energy U, the reactor volume V, the mass of each species, and the coverage of surface species on walls of the reactor.
Equations from Cantera Source
Volume Equation
The equation for reactor volume is:
$ \frac{dV}{dt} = \sum_{walls} \left( K_w A_w \Delta P + F_w(t) \right) $
where $ K_w $ is a wall factor that controls the dependence of wall velocity on pressure gradient across the wall, $ A_w $ is the wall area, $ \Delta P $ is the pressure difference across the wall (so that the wall velocity is pressure-dependent), and $ F_w(t) $ is an arbitrary, user-specified, time-dependent volumetric source term.
Volume Equation Code
The RHS of the volume change equation, above, is assembled by summing over each wall's $ \dot{V} $. This term is computed in Wall::vdot():
/**
* The volume rate of change is given by
* \f[ \dot V = K A (P_{left} - P_{right}) + F(t) \f]
* where \f$ F(t) \f$ is a specified function of time.
*
* This method is used by class Reactor to compute the
* rate of volume change of the reactor.
*/
doublereal Wall::vdot(doublereal t)
{
double rate1 = m_k * m_area *
(m_left->pressure() - m_right->pressure());
if (m_vf) {
rate1 += m_area * m_vf->eval(t);
}
return rate1;
}
The actual computation of these terms, and their assembly into the RHS of the volume change equation above, happens in the Reactor::evalEqs() function:
m_vdot = 0.0;
// compute wall terms
size_t loc = m_nsp+2;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
int lr = 1 - 2*m_lr[i];
double vdot = lr*m_wall[i]->vdot(time);
m_vdot += vdot;
[...]
// volume equation
ydot[1] = m_vdot;
Energy Equation
To track the energy state of the gas in the batch reactor, an internal energy differential equation is solved. This differential equation is of the form:
$ \frac{dU}{dt} = -P \frac{dV}{dt} - \sum_{walls} \left( h_w A_w \Delta T + G_w(t) \right) + \sum_{inlets} \dot{m}_i h_{i,in} + \left( \sum_{outlets} \dot{m}_i \right) h $
where $ h_w $ is the heat transfer coefficient for the wall, $ A_w $ is the wall area, $ \Delta T $ is the temperature difference across the wall, $ G_w(t) $ is an arbitrary, user-specified, time-dependent energy source term, $ \dot{m}_i $ is the mass flowrate of all streams entering into the reactor, $ h_{i,in} $ is the enthalpy of the streams entering into the reactor, and $ h $ is the enthalpy of any stream leaving the reactor.
Energy Equation Code
As with the volume equation, the wall heat flux terms are computed in the Wall class, specifically Wall:Q():
/**
* The heat flux is given by
* \f[ Q = h A (T_{left} - T_{right}) + A G(t) \f]
* where h is the heat transfer coefficient, and
* \f$ G(t) \f$ is a specified function of time.
*/
doublereal Wall::Q(doublereal t)
{
double q1 = (m_area * m_rrth) *
(m_left->temperature() - m_right->temperature());
if (m_emiss > 0.0) {
double tl = m_left->temperature();
double tr = m_right->temperature();
q1 += m_emiss * m_area * StefanBoltz * (tl*tl*tl*tl - tr*tr*tr*tr);
}
if (m_qf) {
q1 += m_area * m_qf->eval(t);
}
return q1;
}
Species Equations (Gas)
The mass of each species $ M_k $ in the batch reactor is tracked with a differential equation of the form:
$ \frac{dM_k}{dt} = MW_k \left[ V \dot{\omega}_k + \sum_{walls} A_w \dot{s}_{k,w} \right] + \sum_{inlets} \dot{m}_i Y_{k,i,in} - \left( \sum_{outlets} \dot{m}_{o} \right) Y_k $
where:
- $ MW_k $ is molecular weight of species k
- $ \dot{\omega}_k $ is the net production rate of species k
- $ \dot{s}_{k,w} $ is the production rate of species k by the wall (surface)
- $ \dot{m}_i $ is mass flowrate of stream i entering reactor
- $ Y_{k,i,in} $ is mass fraction of species k in stream i entering reactor
- $ \dot{m}_{o} $ is mass flowrate out of reactor
- $ Y_k $ is mass fraction of species k in batch reactor (in stream exiting reactor)
The four terms going into the RHS of the species mass balance are:
- Volumetric rate of production
- Surface rate of production
- Mass flow into reactor
- Mass flow out of reactor
The net production rate of species k can be further broken down as:
$ \dot{\omega}_k = \sum_{j=1}^{N_{rxns}} \nu_{jk} r_j $
where $ \nu_{jk} $ is the stoichiometric coefficient of species k in reaction j, and $ r_j $ is the molar production rate for the reaction j.
Species Equation Code
The species equation RHS is assembled, as with the other equations, in Reactor::evalEqs().
The volumetric production rate (first term in the equation above) and the inflow and outflow terms are computed after the loop over each wall:
/* species equations
* Equation is:
* \dot M_k = \hat W_k \dot\omega_k + \dot m_{in} Y_{k,in}
* - \dot m_{out} Y_{k} + A \dot s_k.
*/
const vector_fp& mw = m_thermo->molecularWeights();
if (m_chem) {
m_kin->getNetProductionRates(ydot+2); // "omega dot"
} else {
fill(ydot + 2, ydot + 2 + m_nsp, 0.0);
}
for (size_t n = 0; n < m_nsp; n++) {
ydot[n+2] *= m_vol; // moles/s/m^3 -> moles/s
ydot[n+2] += m_sdot[n];
ydot[n+2] *= mw[n];
}
Here, the line
m_kin->getNetProductionRates(ydot+2); // "omega dot"
gets the rates of production corresponding to each species, following the equation:
$ \dot{\omega}_k = \sum_{j=1}^{N_{rxns}} \nu_{jk} r_j $
The surface production rate (second term in the equation above) is computed in the loop over each wall:
// compute wall terms
size_t loc = m_nsp+2;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
[....]
if (surf && kin) {
[...]
kin->getNetProductionRates(DATA_PTR(m_work));
size_t ns = kin->surfacePhaseIndex();
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
double wallarea = m_wall[i]->area();
for (size_t k = 0; k < m_nsp; k++) {
m_sdot[k] += m_work[k]*wallarea;
}
}
}
It's important to notice here that there are two RHS terms being computed:
The first is happening inside a loop over nk. This is the source term for the surface species, and nk is the number of surface species associated with this particular wall. This is used in the surface coverage equation, which is covered below.
The second is happening inside a loop over m_nsp, which is a class-specific variable (hence the m_ prefix) that stores the total number of species in the gas phase. This second term is the surface rate of production that is going into the species balance equation.
Surface Coverage Equation
See also: Cantera/Surface Coverage
In order to model surface reactions, Cantera treats them similar to the way it treats gas phase reactions. In a gas phase reaction, the reaction rate of a simple reaction like:
$ A + B \rightarrow^{k_{ABC}} C $
depends on the concentrations of A and B. Assuming a first order reaction, and calling this the ABC reaction, the reaction rate looks like:
$ r_{ABC} = k_{ABC} n_{A} n_{B} $
where $ n_{i} $ is the molar concentration of i, with units of moles/volume.
This approach can be extended, by analogy, to surface reactions. If a surface reaction proceeds as (asterisk denotes surface species):
$ A + B* \rightarrow^{k_{surf}} C* $
then the reaction rate looks like:
$ r_{surf} = k_{surf} n_{A} n_{B*} $
The problem, though, is that $ [B*] $ can't be defined in the same way as $ [A] $.
The solution is to define $ [B*] $ as a number of "moles" of a surface species, based on two quantities:
- Site density - the number of sites on the surface available to be "occupied" by any surface species
- Surface coverage - the fraction of all sites on the surface that are occupied by a particular species
This makes the total number of "moles" of surface species B:
$ n_{B*} = \theta_{B} \varrho $
where $ \theta_{B} $ is the surface coverage of species B, and $ \varrho $ is the site density of the surface.
The surface coverage, then, is a "mole fraction" of a species on the surface. The surface coverages also evolve with time, due to changing reaction conditions and production rates, so the surface coverages must be evolved in time using a differential equation, just like the species.
The surface coverage equation for species m is given by:
$ \frac{d \theta_m}{dt} = \frac{ \dot{s}_m \sigma_m }{ \varrho } $
where:
- $ \dot{s}_m $ is the surface production rate of species m (the sum of surface reaction rates involving species m)
- $ \sigma_m $ is the total number of surface sites occupied by species m
- $ \varrho $ is the site density of the surface
Surface Coverage Equation Code
The surface coverage differential equation solution technique is to make the surface coverages part of the solution vector that is being advanced in the Reactor class. The Reactor class already solves a coupled set of ODEs for the volume, energy, and species, so the surface coverages come along for the ride.
This can be seen in the Reactor::evalEqs method, where the wall terms are computed:
// compute wall terms
size_t loc = m_nsp+2;
fill(m_sdot.begin(), m_sdot.end(), 0.0);
for (size_t i = 0; i < m_nwalls; i++) {
int lr = 1 - 2*m_lr[i];
double vdot = lr*m_wall[i]->vdot(time);
m_vdot += vdot;
m_Q += lr*m_wall[i]->Q(time);
Kinetics* kin = m_wall[i]->kinetics(m_lr[i]);
SurfPhase* surf = m_wall[i]->surface(m_lr[i]);
if (surf && kin) {
double rs0 = 1.0/surf->siteDensity();
size_t nk = surf->nSpecies();
double sum = 0.0;
surf->setTemperature(m_state[0]);
m_wall[i]->syncCoverages(m_lr[i]);
kin->getNetProductionRates(DATA_PTR(m_work));
size_t ns = kin->surfacePhaseIndex();
size_t surfloc = kin->kineticsSpeciesIndex(0,ns);
for (size_t k = 1; k < nk; k++) {
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
sum -= ydot[loc + k];
}
ydot[loc] = sum;
loc += nk;
[...]
}
}
The for loop over nk, where the terms are being put into ydot, is where the actual RHS for the surface coverage differential equation is constructed. The m_work array is populated with the surface production rates (just a few lines above), and surf->size(k) gets the total number of surface sites occupied by species k, so the expression:
ydot[loc + k] = m_work[surfloc+k]*rs0*surf->size(k);
is constructing the RHS of the equation:
$ \frac{d \theta_m}{dt} = \frac{ \dot{s}_m \sigma_m }{ \varrho } $
Equations from Froment and Bischoff
1D Model
From p. 403. Conservation eqn for steady state, single rxn carried out in cylindrical tube:
1D Material Balance
$ -\dfrac{d}{dz} \left( u_s C_A \right) = r_A \rho_B $
1D Energy Balance
$ u_s \rho_g c_p \dfrac{dT}{dz} = \left( - \Delta H \right) r_A \rho_B - r \frac{U}{d_t} ( T - T_r ) $
1D Momentum Balance
$ - \frac{dp_t }{dz} = f \frac{ \rho_g u_s^2 }{d_p} $
Friction factor dependent on shape factors, &c.
Reynolds number, a, b, hollow rings, packed bed of spheres, etc.
2D Equations
The 2 dimensional equations, accounting for radial mixing, are treated with an effective diffusivity method. The equations become coupled two-point boundary value problems. The governing equations become:
Material Balance
$ \varepsilon D_{ea} \dfrac{d^2 C_A}{dz^2} - u_s \frac{dC_A}{dz} - r_A \rho_B = 0 $
Energy Equation
$ \lambda_{ea} \frac{d^2 T}{dz^2} - \rho_g u_s c_p \frac{dT}{dz} + \left( - \Delta H \right) r_A \rho_B - \frac{4U}{d_t} \left( T - T_r \right) = 0 $
Boundary Conditions
$ u_s ( C_{A0} - C_A ) = - \varepsilon D_{ea} \frac{dC_A}{dz} \qquad z = 0 $
$ \rho_g u_s c_p (T_0 - T) = - \lambda_{ea} \frac{dT}{dz} $
$ \frac{d C_A}{dz} = \frac{dT}{dz} = 0 \qquad z=L $