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.

Also, see this note about how to set up of computationally efficient PDE coefficients.

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.

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).

Related Tech Notes

▪
• Mass Transport Analysis
▪
• Mass Transport Model Collection
▪
• PDEModels Overview