Microscale Simulation of Catalyst Deactivation
| Introduction | Solve the PDE Model |
| Mass Transport Model | Appendix |
| Domain | Nomenclature |
| Material Parameters | References |
| Initial and Boundary Conditions |
Introduction
A catalyst is a substance that increases the rate of a chemical reaction without being consumed itself. Although it is not consumed the catalyst can be deactivated by secondary processes. Two types of catalyst deactivation are common. First, the loss of a catalyst's reactivity, the ability to increase the rate of a reaction. Second, the loss of selectivity, the ability to direct a reaction to yield a particular product over time [1]. To avoid insufficient catalyst activity and selectivity the deactivation of catalysts should be understood and monitored. This study is to simulate a gas-solid catalyst reaction within a porous structure and model the deactivation process over time.
The type of catalyst deactivation we look at in this study is due to coke accumulation [2]. In this case coke accumulates over catalyst sites, 
, and reduces the catalyst's reactivity. In this model a gas-phase reactant, 
, is passed over a porous structure of a catalyst particle, and is converted into a gas-phase product, 
. This product will then react with the catalyst as the secondary process, and form coke, 
, on the catalyst sites. This coke accumulation eventually deactivates catalyst sites and reduces the rate of the reaction possible at the sites.
is the reactant,
is the catalyst sites,
is the product,
is the accumulated coke,
are the reaction rate constants.
To build a realistic geometry that represents the porous feature of a catalyst microstructure, we will use a microscopy image of a catalyst particle [3].
The catalyst particle is roughly assumed to be spherical. As such the reactant's concentration will decrease with radius 
. Because of symmetry it is then sufficient to use a 2D cutout of the catalyst particle for the simulation.
In this model we assume that there is no external fluid flow involved, which means the reactant 
is transported by diffusion only; there is no convective mass transfer.
To measure the level of the catalyst deactivation, we can compute the amount of coke 
that accumulated on the catalyst sites 
. The net increase in the coke concentration is computed by integrating over the porous structure of the catalyst particle.
Here 
denotes the local concentration of coke 
and 
denotes the mean concentration of coke 
.
The symbols and corresponding units used here are summarized in the Nomenclature section.
Please refer to the information provided in the "Mass Transport" theory manual for a more general theoretical background for mass transport/chemical reaction analysis.
Mass Transport Model
The non-conservative mass transport equation (4) is used to solve for the concentration distribution 
in a mass transport model:
is the concentration of the transported species 
,
is the species diffusivity 
,
is the rate of mass production/consumption, the volumetric reacting rate of the species 
.
The kinetic model of the catalyst reaction we wish to model is given by:
Here 
and 
are the reaction rate constants of this two-steps process. First we will simulate the concentration evolution of reactant 
, catalyst sites 
and the resulting product 
. Once we have computed the concentration profile of 
, 
and 
, we can proceed to solve for the accumulation of coke 
based on the second step of the reaction in (5).
During the first step the concentration 
of the reactant 
denoted as 
will reduce while the product concentration 
increase. It is important to note that the concentration of catalyst sites 
will only get consumed in the second reaction, resulting in a raise of the coke concentration 
.
We assume the catalyst reaction to be a first-order reaction [6], that is, the reaction rate 
of each reacting species will be linearly proportional to the concentration of the involving reactants:
The values of rate constants 
and 
were estimated from the literature [7].
In this model we assume that there is no external fluid flow involved, which means the mass convection term: 
in (8) is set to zero and the mass is transported by diffusion only.
The reaction terms could very well also me modeled with a "MassReactionRate" model parameter in place of the "MassSource" term. Then a dependent variable needs to be factored out of the reaction rates and the sign of the equations negated. The choice is purely a matter of taste.
Domain
In this model we assume that the catalyst reaction occurs homogeneously from the particle surface into the core of the particle. As such the variation of the species concentration in the azimuthal direction 
and the polar direction 
can be ignored. Then a two dimensional model is sufficient to represent the 3D catalyst particle. Here we make use of the lower-left quarter of the catalyst particle as the modeling domain 
.
To build a realistic geometry that represents the catalyst microstructure, we will use a microscopy image of the catalyst [9]. Note that the image needs to be binarized in order to distinguish the gas phase region (shown in white) and the solid phase region (shown in black).
The binarized image has a finely resolved catalyst boundary. When we discretize this image into a mesh a large number of elements is expected. Since a large number of elements will result in a longer run time and a larger memory usage one may consider to defeatured the domain. Defeaturing is done by making use of a less accurate catalyst boundary while still preserving the important characteristics of the original structure. The details about defeaturing the domain and its construction is presented in the appendix section: Construction of De-featured Modeling Domain. For the sake of accuracy we make use of the original region.
Next, we convert the binarized image into a mesh domain in order to apply the finite element method.
To model the catalyst reaction in the realistic scale, the domain is resized based on the radius of the catalyst particle measured in [10]: 
.
To generate the full element mesh, we have to set up element markers to differentiate the gas phase region and the solid phase region. These markers can then be used to specify, for example, the species diffusivity
in different regions of the simulation domain.
For markers to be attributed to different subregions, coordinates that lie in those subregions need to be specified. Here we extract the coordinates of the solid phase region (catalyst sites 
).
Here the gas phase domain 
is shown in white and the solid phase domain 
is shown in blue. The two subregions are separated by internal mesh boundaries.
More information on markers and their generation in meshes can be found in the Element Mesh Generation: Markers documentation.
Material Parameters
In this model the diffusivity for the reactant 
and the product 
are different between the gas and solid domains. The diffusivity of the catalyst sites 
is set to zero to represent its stationary feature.
Initial and Boundary Conditions
In this model we are simulating the catalyst deactivating reaction for five seconds.
At 
the initial concentration for the reactant 
and the product 
are zero throughout the domain 
. The catalyst sites 
is confined in the solid phase 
with an initial concentration of 
.
A simple way to model the discontinuous initial concentration 
is by using the If statement with element markers:
However, for the default second order mesh a discontinuous initial condition may result in overshoot/undershoot values for 
on the gas-solid interface due to the quadratic interpolation. This issue is explained in a separate tutorial: Finite Element Method Usage Tips.
To improve the accuracy we can use a linear interpolation instead. To do so we create a first order mesh and extract the original (second order) interpolation result on first order nodes. Then we compute the linear interpolation based on these values.
There are three types of the boundary conditions involved in this mass transport model: the concentration boundary condition, the surrounding flux boundary condition, and the symmetric boundary condition.
The following illustration shows the locations where these boundary conditions are applied.
A default symmetric boundary condition is implicitly applied on the axis of symmetry.
A concentration boundary condition is used to model the supply of the reactant 
on the open surface. To avoid a discontinuous jump of reactant at the boundary a helper function that smoothly ramps up the reactant supplied over time is created.
During the catalyst reaction the product 
is free to diffuse out of the simulation domain through the open surface. This is modeled by a surrounding flux boundary condition.
and the mass transfer coefficient 
on the open surface for the surrounding flux boundary condition for the product 
. Solve the PDE Model
Concentration Field for G, S and P
First we solve for the concentration field for the reactant 
, the catalyst sites 
and the product 
.
As mentioned in the previous section, due to the large number of elements used we expect a longer run time and a larger memory usage to find a solution than in typical application examples shown. The number of the degrees of freedom is proportional to the number of elements and gives a measure of the effort needed to solve a PDE.
Since we have 3 coupled PDEs are solving for 
(reactant 
, product 
and catalyst sites 
) nodal values at each time step.
We solve for almost 2 million degrees of freedom.
During the reaction the reactant 
gradually diffuses into the porous structure of the catalyst particle, and is converted into the product 
. This product 
then reacted within the porous structure and forms the coke 
on the catalyst sites 
. The coke accumulation eventually deactivated catalyst sites by decreasing the concentration of the catalyst sites.
The plots above are generated with a call to Rasterize. This is to reduce the disk space this notebook requires. To obtain high quality graphics comment out the call to Rasterize.
Concentration Field for B
Once we have computed the concentration field of the product 
and the catalyst 
, we can proceed to solve for the accumulation of coke 
based on the reaction equation:
The kinetic model for the accumulated coke 
is given by:
To represent the stationary feature of the solid-phase coke 
, the diffusion coefficient is set to 
and the equation (11) simplifies to:
The coke 
gradually accumulates throughout the porous structure of the catalyst particle over time. To measure the level of the coke accumulation, the increase in the mean coke concentration: 
is computed. This is done by integrating over the solid phase domain 
:
Due to the catalyst reaction, the coke 
has accumulated to 
from 
to 
.
Appendix
Construction of Defeatured Modeling Domain
Due the large number of elements 
used, it took a significant amount of time/memory to solve this PDE model. To reduce the number of elements and the time/memory usage, we can defeatured the domain instead, that is, we model a less accurate catalyst boundary but still preserving the important characteristics of the original structure.
To build a defeatured domain we will smooth out the catalyst boundary in the binarized image with DeleteSmallComponents.
Here 
and 
are parameters that control the level of defeaturing, and should be carefully chosen to preserve the accuracy of the modeling domain.
Note that the defeatured domain has a smoother catalyst boundary, which will significantly reduced the number of elements during the domain discretization.
Note that the mesh elements has been greatly reduced from 322k to 100k, leading to a lower time/memory consumption when solving the model.
The modeling result of the defeatured domain is given by:
Compared with the result from the original domain: here and here, it is seen that the catalyst reaction took a slightly longer time to propagate to the inner region through the porous structure. This is because some minor channels within particle has been smoothen out and treated as a solid phase region with much smaller diffusivity, which makes it more difficult for the reactant 
to pass through.
Nomenclature
References
1. Bartholomew, Calvin H. Mechanisms of Catalyst Deactivation, Applied Catalysis A: General 212 17–60 (2001).
2. Wolf, E. E. and Alfani, F. Catalysts Deactivation by Coking, Catalysis Reviews: vol 24 329-371 (1982).
3. Ciesielski, P. N., Robichaud, D., Donohoe, B., Nimlos M., Microscale Simulation of Catalyst Deactivation during Gas-Phase Upgrading of Biomass Pyrolysis Vapors, Biosciences Center, National Renewable Energy Laboratory (2015).