Single catalyst particle growth modeling in thermocatalytic decomposition of methane

ThermoCatalytic Decomposition of methane (TCD) is studied as a method to convert natural gas into hydrogen and functional carbon. In these processes the carbon typically formed on top of a catalyst phase leading to particle growth. Therefore, the development of a particle growth model is necessary to understand the limitations of thermocatalytic decomposition of methane and to assess optimal parameters and process conditions. In this paper, a particle growth model is presented to describe the growth of functional carbon on the catalyst particle. This coupled model requires kinetic equations and information on deactivation rates which have been studied from literature. The morphology of the particle changes due to carbon formation, which leads to eventual deactivation. Therefore, these kinetic expressions are coupled to a particle growth model based on the analogy with the growth of particles in polyolefin production. To combine the effects of particle growth, kinetics, and internal heat and mass transfer, the Multi-Grain Model (MGM) was used. Results confirm that with the currently available catalysts the carbon yield is not affected by heat and mass transfer limitations, however, with the availability of more active catalysts these limitations will become important. Temperature, however, has a significant role in that it regulates the kinetic rate and thus growth rate, which in turn influences the catalyst deactivation. The optimum temperature for the production of nano-carbon, within a reasonable process time, therefore sensitively depends on the choice of catalyst.


Introduction
Hydrogen can be produced through different processes from different feedstocks, such as steam methane reforming, water splitting, and thermocatalytic decomposition of methane. Steam methane reforming coupled with CO 2 capture and storage (CCS) technologies are the most known and investigated methods in the field of low carbon footprint hydrogen production. However, the separation of produced CO 2 and handling and storage leads to costs for gas separation and storage management. Water splitting is an energy and capital intensive process and increases the final price of hydrogen. By contrast, methane decomposition to functional carbon materials and hydrogen has advantages to the alternative processes such as the elimination of additional purification/separation units and production of valuable carbon nanomaterials (tubes or fibers) instead of CO 2 . The potential applications of carbon nanomaterials in semiconductors, additives to building materials, energy storage, and catalytic materials due to their unique physicochemical properties such as high conductivity, high tenacity and mechanical strength, high specific surface temperature range and is also more likely to deactivate. However, the addition of a second (or even a third) metal such as iron or copper to the nickel, has shown improved stability at higher temperatures and less deactivation [11][12][13][14].
Despite the high potential of TCD for producing carbon nanomaterials and CO 2 -free hydrogen, it is greatly restricted for industrial applications due to inadequate productivity, uncertainties of process Density of the grain (kg∕m 3 ) Density of produced carbon layer, including porosity (kg∕m 3 ) performance and operational challenges coming from carbon formation [2]. Therefore, in addition to studies on catalysts, the TCD reactor and process has to be developed, designed and controlled thoroughly to become feasible at industrial scale. For a rational reactor and process design, modeling and experimental studies can provide the required understanding and basic data for this. This understanding facilitates identification of optimal process conditions for maximum carbon nanomaterial production. Modeling of the catalyst performance as function of equivalent process time is critical for understanding and predicting the product and catalyst evolution in the reactor. This performance can be expressed with the ratio of the mass of produced carbon to the mass of fresh catalyst used, called carbon yield (Eq. (2)) and the change in catalyst particle size and density.
In the literature, including the work of Ashik some studies have been reported on modeling at the molecular scale [2,15,16], which helps scientists in catalyst evaluation and to develop a microscopic level of understanding of the complete chemical transformation. Although these models are helpful in the understanding of reaction mechanisms, a model that properly describes the behavior of a catalyst at the macro level has not been reported to the best of our knowledge. In particular, the formation of a functional carbon layer onto the catalyst phase leads to particle growth [2]. The particle size and thus growth is a crucial parameter in designing a reactor. In the present work, for the first time, the multi-grain model (MGM) based on the analogy with the growth of polyolefin particles is developed to describe the macroscopic behavior of growing particles in TCD of methane. The model couples different phenomena involved inside the catalyst particle (which is called macroparticle in this study) such as heat and mass transfer and chemical reaction. A short review of kinetic studies in the literature is provided in Section 2. In Section 3, the model description is presented. The model validation and its reactor predictions assessed in Sections 4 and 5 respectively.

Short review of kinetics in literature
In TCD, the activity of the catalyst and kinetic rate of reaction decrease over time due to deactivation. The actual rate of reaction at time > 0 is described using two parameters: the initial reaction rate and a time-dependent deactivation factor.

Initial reaction rate
The initial reaction rate and the reaction mechanism have been studied extensively [2]. Douven et al. and Yadav et al. proposed reaction rate equations for carbon nano-tubes (CNT) production by TCD of the methane which is only dependent on methane concentration [3,17], (see Eq. (3) below). Yadav found out that multi-walled CNT is produced with a different kinetic rate than single-walled CNT, both only depend on the concentration of methane [18], (see Eq. (4)).  proposed Eq. (5) for the initial reaction rate which also does not depend on the hydrogen concentration [19].
The majority of studies have revealed that the hydrogen concentration has a negative effect on the initial reaction rate due to thermodynamic factors equilibrium and occurrence of the reverse reaction [5,7,[20][21][22][23]. The kinetic models presented in Eqs. (3)-(5) must therefore be regarded as a simplified form of the actual kinetics. The equations Eqs. (6)-(8) that involve the effect of the hydrogen concentration have a very similar form, but differ mostly in the expression in the denominator. Amin et al. [7] and Snoeck et al. [21] derived Eq. (6) while Borghei et al. [22] suggested different powers for H 2 and CH 4 , Eq. (7), which may be due to the use of a different catalyst and the specific operating conditions used in their studies [7,21,22]. Saraswat et al. [5] have reported a more extended form of the reaction rate, Eq. (8), which includes effects of hydrogen and methane partial pressures in Langmuir-Hinshelwood type of expressions [5].
According to the literature [5,7,23] The most likely mechanism of the reaction is based on molecular adsorption of methane as the first step, followed by a series of dehydrogenation reactions that are taking place one by one until it ends with separate adsorbed carbon and hydrogen atoms. Detachment of the first from 4 is known to be the slowest and the rate determining step. Every two adsorbed hydrogen atoms form a single 2 molecule, which is released from the catalytic surface. The carbon atom, however, can diffuse into the nickel catalyst and either forms nanomaterials or encapsulates the active site.

Deactivation factor
The ratio of reaction rate at time to the initial reaction rate (Eq. (9)) is called the deactivation factor. The deactivation factor expresses the stability of the catalyst over time, which is a crucial factor in order to obtain a high carbon yield.
Several different empirical or semi-empirical equations for the deactivation factor a have been defined. Borghei [22] proposed Eq. (10) for the deactivation factor, where b, c and d are constants and is defined by an Arrhenius type of temperature dependency. Douven [3] reported that deactivation is reversible and probably due to the formation of amorphous carbon and encapsulation of active sites. Douven used a sigmoid Eq. (11) as the deactivation factor, where parameter b is assumed to have only temperature dependency however parameters 0 and c decrease slightly as the methane partial pressure increases.
Eq. (12) is proposed by Amin et al. [7] and is based on the proposed mechanism in 2.1, mass balance of species on the surface of the catalyst and the assumption that all the reaction steps except one are in equilibrium. So, the concentrations of intermediate species are negligible. All parameters ( , , , , 4 and , 2 ) are temperature dependent following the Arrhenius equation and are determined by fitting the expression to experimental data.

Model description
For our model of the TCD reactor and the particle growth over time, the following physical phenomena have to be taken into account: 1. The transport of species into and out of the particle, being diffusion of methane into the catalyst pores and diffusion of hydrogen out of the catalyst pores. 2. The heat transfer into the macroparticle to provide the heat for the strongly endothermic reaction. This includes the transfer from the bulk gas to the external surface of the macroparticle, conduction within the macro macroparticle as well as conduction within the grains. 3. The decomposition of methane on the active sites of the macroparticles into hydrogen and solid carbon. 4. The accumulation of solid carbon onto the catalyst, with consequential growth of the catalyst and deactivation of the catalyst.
These phenomena are very similar to the olefin polymerization, which also experiences particle growth of the macroparticle due to solid product formation. For the solids formation and particle growth different modeling approaches have been developed [24,25]. For our TCD process the most appropriate model, capable of modeling particle growth and convenient for further development to include fragmentation [26], is the Multi-grain model (MGM).
MGM is based on two assumptions. Firstly, the macroparticle is spherical with only profiles in the radial direction. So, it is a 1D model in the radial direction. Secondly, the macroparticle is composed of layers of identical non-porous grains (microparticles) with active sites on their surface [25][26][27][28][29]. The growth and evolution of the macroparticle are due to the accumulation of produced carbon on the surface of grains. Fig. 1 shows the schematic of the both concepts of macroparticle and internal grain layers before and after the reaction takes place. Fig. 1(a) shows the schematic of a fresh porous macroparticle which consists of layers of non-porous grains that are illustrated as black spheres. Fig. 1(b) shows a circular sector of the same macroparticle after entering the reactor. The gray shell around each grain is the produced carbon on the grain.

Mass transfer
Methane and hydrogen diffusion through the pores of the macroparticle and through the layer of accumulated carbon surrounding the catalyst fragments is modeled at two different scales. At the scale of the macroparticle the process is modeled by the diffusion-reaction equation in spherical co-ordinates: Where is the molar concentration of a component, is the radial position in macroparticle and is the effective diffusivity of the considered component. ( , ) is the average rate of reaction at a given radial position: Where is the porosity of catalyst, , is the number of grains at radial position and 0 is the radius of the core of the grains. Eq. (13) can be solved with the following boundary and initial conditions: is the external mass transfer coefficient of the macroparticle, is the external concentration and 0 is the initial concentration in the particle.
The concentration at the grain scale is also modeled by the diffusion equation in spherical coordinates, Eq. (18). Considering the assumption that the core of the grains are non-porous, so there is no hydrogen or methane inside the core of the grains. Therefore, the boundary conditions, Eqs. (19) and (20) are defined at the surface of the core and outer surface of the grain particle.
Where is the radial position in the grain, is the diffusivity in the product layer around the grains, is the radius of whole-grain and , is the position of grain in the macroparticle (Fig. 1) and is the unit conversion factor. The initial condition is the same for the macroparticle (Eq. (17)). Eqs. (13)- (17) are solved only once during each time step for the macro macroparticle, while Eq. (18) and its initial and boundary conditions are solved for all the internal grain layers at each time step.

Heat transfer
The heat transfer mechanism at both scales is based on conduction. The temperature profile in the macro-particle is calculated by Eq. (21): Where ℎ , and are the heat conductivity, the density and the specific heat capacity of the macroparticle respectively, is the heat of reaction and ( , ) is obtained from Eq. (14). Eq. (21) can be solved with the following initial and boundary conditions (Eq. (22) to (24)): ℎ is the external convective heat transfer coefficient outside the macroparticle which can be estimated from the Gunn correlation, 0 and are respectively the initial temperature and the temperature outside the macroparticle.
The radial temperature profile in the grains can be obtained from the heat conductivity equation in spherical co-ordinates for the whole domain of the core of the grain and the layer of carbon product. However, the heat conductivity, density and specific heat capacity ( , and ) have different values in the core of the grains and in the product layer.
One can notice the analogy between mass and heat transfer at both scales of the macroparticle and the grains.

Reaction kinetics
In this study, Eq. (6) and (12) and corresponding parameter values are provided in Table 1 and are used in the model as the initial reaction rate and deactivation factor of reaction [7]. Carbon formation increases the radius of the grains and consequently the size of the macroparticles. The growth rate of the radius of the grains is calculated by Eq. (28) and the growth rate of the macroparticle equals the summation of the growth rate of all grain layers, Eq. (29).

Verification of the model
The performance and the results of the model are validated by comparing its results with analytical solutions and results obtained from an independent PDE solver. Two limiting cases are used to verify the implementation of the model. In the simplified case 1, it is assumed that there is only one layer of micro grains, the mass transfer limitation is low, and the reaction is first order without deactivation. In this case, the mass of carbon produced is calculated by Eq. (30): Where 0 is the radius of the core of micro grains, is the number of micro grains in the only available layer in the macroparticle, is the time passed since the start of the reaction, is the kinetic coefficient of the first-order reaction per surface area of grain core (mol∕m 2 ) and 4 is the molar mass of methane. Fig. 2 illustrates the high accuracy of MGM in case 1, by showing that the MGM results matches Eq. (30).
In case 2, again it is assumed that the reaction is first order in methane and independent of the hydrogen concentration without any deactivation ( = . 4 ). The second assumption is that the reaction takes place uniformly in the macroparticle. Finally, it is assumed that the physical properties of the macroparticle do not change with time as the reaction proceeds. In these conditions the methane concentration profile inside the particle can be calculated at any time by solving Eq. (31) and its associated initial and boundary conditions.
( , = 0) = 0 ( = 0, ) = 0 The set of equations can be rewritten in dimensionless from via definition of the following dimensionless quantities: is a Thiele modulus and represents the ratio of reaction rate to diffusion rate. By using these dimensionless quantities in Eqs. (31), they change to: Eq. (36) has been solved with an independent PDE-solver (Matlab pdepe solver). Fig. 3 compares the results of MGM with the Matlab solver. As is clearly illustrated, the results of the MGM are in almost perfect agreement with the Matlab solver over the time, from start of the reaction till the steady state condition has been reached (which is about one minute in this case). This confirms that the model works correctly and there are no errors in the calculation methods.

Results and discussion
The results of MGM simulations are evaluated by means of a sensitivity analysis and the assessing importance of different parameters and operating conditions in TCD. The temperature range used in our models is limited to the operating conditions that are used to derive the kinetic coefficients of Eqs. (6) and (12) (temperature: 500-650 • C and maximum hydrogen fraction: 10%) to ensure of the validity of results [7]. Table 2 provides the most important parameters used in the model for the following cases, unless otherwise stated or the importance of the parameter is evaluated.

Mass and heat transfer limitations
To assess the impact of internal mass transfer limitations, the diffusivity was altered to one thousand times higher and lower values than the (base case) values stated in Table 2. The results are summarized in Fig. 4 and reveal that lowering the internal mass transfer rate does not affect the carbon yield. Hence, for the conditions of Table 2 the effect of diffusion limitation is negligible. On the other hand, if the mass transfer limitation increases to one thousand times higher, the carbon yield decreases by about 35%.   The importance of internal diffusional resistance is also confirmed by the Weisz-Prater criterion (Eq. (37)) that estimates the importance of the diffusion on the reaction rates in heterogeneous catalytic reactions [30]. In the normal case, = 9.02 × 10 −4 ≪ 1 which means that internal mass transfer does not influence the production rate of carbon. However, in the case with one thousand times lower effective diffusion coefficient, = 0.902 which implies that for such low diffusion coefficients, internal mass transfer limitation is not negligible anymore compared to the reaction rate. The same procedure was applied to assess the role of the internal heat transfer limitation, by changing the thermal conductivity of the solid material composing the macroparticle. The results are presented in Fig. 5. Since the carbon yield is not affected by lowering the heat transfer limitation, for the conditions of Table 2 the heat transport resistance is negligible in comparison with other factors. However, in the case with 1000 times higher heat transfer limitation, the carbon yield is decreased dramatically. In this case, the temperature in the macroparticle increases relatively slowly, and as a result the reaction rate and therefore the slope of the curve increases gradually.
In addition, there can prevail external heat and mass transfer limitations in the thin film around the macroparticle. However, it was observed that even with the highest external heat and mass transfer limitation, meaning a macroparticle in a stagnant gas phase ( = 2 and ℎ = 2) and with larger particle sizes (1000 μm), the production rate of carbon is not reduced.
These observations regarding mass transfer importance and their effect on the TCD process are in agreement with literature findings derived from both experiments and the Weisz-Prater criterion [5,22]. Fig. 6 illustrates in logarithmic scale how much the carbon yield changes if the initial reaction rate changes by a factor 1000. Reduction of the reaction rate leads to decrease by a factor 1000 in carbon yield. This finding is another confirmation of the fact that the reaction is the rate-determining step compared to the mass and heat transfer limitations (Section 5.1). On the other hand, if the reaction is one thousand times faster, the carbon yield increases around 450 times. This means that in this case mass and heat transfer limitations become also important which again is in agreement with Section 5.1. Fig. 6. The effect of initial reaction rate when it changes one thousand times. In the case with one thousand times faster kinetics, the carbon yield is increased about 450 times, while the change for the case with one thousand times slower kinetics the difference is about one thousand times.

Bulk gas concentration effect
As can be seen in Fig. 7 adding inert gas (which means lowering the methane fraction) decreases the carbon yield. On the other hand, for a given fraction of methane, increasing the fraction of hydrogen leads to lower initial reaction rate and higher durability of the catalyst against deactivation. These effects are presented in Fig. 8. These two effects together lead to higher carbon yield, however, in comparison a pure methane feed yields a higher amount of carbon in a shorter amount of time.

Temperature effect
Temperature has two opposing effects in the TCD process. On one hand, higher temperature leads to a higher initial reaction rate and therefore a higher carbon production rate. On the other hand, increasing the temperature results in faster deactivation of the catalyst and lowers the final carbon yield. These two phenomena are illustrated in Fig. 9. At high temperatures deactivation proceeds more suddenly instead of gradual deactivation at lower temperatures. Thus, curves of 600 • C and 650 • C have a very short flat part, rather than longer flat tail.
Increasing the temperature between 500 • C to 650 • C leads to a lower carbon yield due to the increased deactivation rate. However, it should be noted that this slightly lower amount of carbon is produced in a significantly shorter period of time. Therefore, in the examined conditions and with the used kinetic model, the optimum operating conditions will depend on economic considerations. It should be noted that using a different catalyst (and as a result, different kinetic models) may change this optimum condition.

Effect of number of micro grain layers
The number of micro grain layers in the macroparticle is a model parameter that is not straightforward to measure or estimate, as previous parameters were. As Fig. 10 shows, this number has a significant impact on the carbon yield. The effect of the number of grain layers is not linear and becomes stronger with an increase in the number of layers. Although physically the number of micro grain layers can be translated to the specific surface area of the macroparticle, the internal structure of an actual catalyst particle is more complex than the structure defined by many layers of identical spheres. Therefore, the number of grain layers will be used as the tuning parameter of the model against validated data.

Conclusions
A Multi-Grain Model has been developed to model the heat and mass transfer inside macroparticles coupled with the decomposition reaction of methane. The reaction rate model and deactivation factor from Amin [7] are used, however, the model is suitable for the use of other kinetic models which can be easily accommodated.
The effect of operating conditions and model parameters has been investigated by sensitivity analyses and it was found that the heat and mass transfer rates do not limit the carbon production rate and consequently the reaction is the rate-limiting step of the process. This fact is in agreement with experimental findings respected in the literature. However, if a catalyst is made with one thousand times higher ratio of kinetics rate to the mass and heat transfer rates (either by increasing Fig. 9. The effect of operating temperature on the carbon yield. Higher temperature leads to higher initial reaction rate (initial slope of the curves) and faster deactivation (end point of curves). the reaction rate or decreasing the mass and heat transfer rates), the heat and mass transfer limitations will affect the final carbon yield.
It was found that, the presence of hydrogen causes a decrease in the reaction rate, however a higher carbon yield is achieved due to delayed deactivation of the catalyst. Moreover, increasing the operating temperature leads to a faster initial reaction rate and faster catalyst deactivation and hence an optimal, process dependent, temperature exists.
In the future, it would be interesting to conduct experimental tests to tune, validate and further develop the model. The findings of the current article can be used in CFD models and enable researchers and industry to model, design and predict the behavior of multiple particles in the fixed or fluidized bed reactors employing TCD.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.