Mathematical modelling of the carbon nanotubes synthesis by methane pyrolysis on a copper catalyst based on the approximation approach

. We present the mathematical model concepts for the synthesis of carbon nanotubes by catalytic pyrolysis of methane, based on the analysis and approximation of the dependence of the process stages on the gas phase composition and temperature. The process was considered to consist of three successive stages that include the linear stage of initial growth with a constant specific rate, the stage of growth retardation due to the catalyst active site deactivation, and the stage close to linear, occurring after equilibrium of the accompanying sorption-desorption processes on the active sites is reached. Based on the available experimental studies on a specific catalyst, a step-by-step search for mathematical patterns describing the selected phases of a typical nanotubes growth curve was performed. We have determined and confirmed the criteria for the dependence of the proposed equations parameters on the gas phase composition. The article also discusses the advantages of the proposed approach for modeling complex processes in comparison with traditional physicochemical approaches.


Introduction
Carbon nanotubes (CNTs) are of great interest because they provide potential for creating innovative products in many areas [1][2][3][4]. By directly controlling the reactions and interactions between atoms, it is possible to produce various nanomaterials and control their properties. CNTs have a unique combination of mechanical, thermal, electrical, chemical and biological properties, which can be varied by changing the synthesis conditions, e.g. the length, diameter, chirality, number of layers, and other morphological features of nanotubes depend on the synthesis conditions [5][6][7][8]. There are several main methods for the synthesis of CNTs, but the method of chemical vapor deposition (CVD) is one of the most common and promising [9][10][11][12]. Using mathematical modeling approaches, we can describe the synthesis of CNTs with a detailed analysis of the process and determine the link between the variable process parameters and the mass yield of CNTs [13][14][15]. The parameters affecting the yield of nanotubes include the ratio of gas components of the feed stream, the catalyst metal elements, and the temperature in the reaction zone [16]. However, an incomplete set of experimental data often serves as a limitation for complex analysis.
The main purpose of this study is to analyze the kinetic curves of CNTs growth and mathematically determine the impact of temperature, methane and hydrogen concentrations on kinetic curves.

Materials and methods
The initial data for the development of a mathematical model were taken from the results of experimental studies on the CNTs synthesis by the CVD method, presented in [16]. The experiments were carried out on the Co-Cu/CDC catalyst (CDC -cellulose-derived carbon). The nominal metal content in the catalyst was Co0.05Cu0.0135CDC0.9365 with respect to the initial weight of cellulose. The gas phase consisted of methane (CH4), hydrogen (H2), and nitrogen (N2). The experiment was performed in a tubular reactor, the catalyst was placed on a suspended substrate, which was inert to the catalyst, then the substrate with the catalyst was lowered into the reaction zone heated to the required temperature. The synthesis began from the moment the gas flow was supplied to the reaction zone. The gas flow rate was controlled with rotameters. Detection of nanotube growth was carried out with torsion balance holding the substrate and catalyst on nichrome filaments. After the completion of the experiment, the reactor was cooled with a stream of argon. The required concentrations of CH4 and H2 were set by mixing the volumetric flows of gases. Figure 1 shows a diagram of the gas phase compositions specified for the experiments, presented in [16]. The study describing influence of the gas phase on the kinetics of CNTs growth was performed at a constant temperature of 750 ºC due to the high carbon productivity and low catalyst deactivation and with a change in the concentration of CH4 and H2 from 7.1 to 42.9 vol%. A horizontal line shows the dependence of nanotube growth on the hydrogen concentration at a constant methane concentration. A vertical line represents the dependence of nanotube growth on the methane concentration at a constant hydrogen concentration. To evaluate the effect of the reaction temperature on the yield and morphology of the grown CNTs, the Co-Cu/CDC catalyst was operated at temperatures from 650 ºC to 975 ºC with a gas phase composition of CH4:H2:N2 = 28.6:14.3:57.1 (vol%). The analysis of the experimental results demonstrates both productivity and carbon growth rate improve with increasing temperature and partial pressure of CH4. However, as the temperature increased above 800 ºC, the carbon yield reached a plateau after 120 min of reaction, that could be associated with partial deactivation of the catalyst either due to carbon encapsulation or steric hindrance. An increase in the concentration of methane leads to a rise in the amount of carbon deposited on the metal-carrier surface due to the carburization of metal nanoparticles and an increase in the number of carbon atoms dissolved in them. In addition, the hydrogen concentration has a significant effect on the growth rate of carbon nanomaterials -the less hydrogen concentration, the faster process. This is due to the adsorption of hydrogen molecules on the active sites of the catalyst and therefore the difficulty of adsorption of methane molecules on the same active sites.
The authors from [16] have also proposed their own mathematical model to describe the process under study, but they have faced difficulties to identify and justify the characteristic trends in the change of their model parameters from the parameters that specify the synthesis conditions. It made the presented model unsuitable for prognosis and capable of describing only partial curves within the framework of the experiments performed. Having analyzed the experience, advantages and disadvantages of the model from [16], the authors of this current article have proposed their own concept for the mathematical description of such processes. The mathematical model is based on the analysis of the impact of temperature and composition of the gas phase on the nanotube growth curve in general and on individual sections in particular. Also, the subsequent search for mathematical dependencies reflecting this effect was performed. On a typical CNTs growth curve of catalytic pyrolysis, three successive stages can be distinguished ( Fig. 2): 1) the stage of initial growth, which can be considered linear with a high degree of accuracy, 2) the stage of catalyst active site deactivation, which leads to the significant slowdown in the growth rate and, as a consequence, a change in the slope of the curve; 3) the stage that occurs after the equilibrium of sorption-desorption processes in the active centers of the catalyst is established. Moreover, stage 3, like stage 1, can be considered linear. It is noteworthy that in a number of experiments, stages 2 and 3 may be absent, which is probably due to the rapid restoration of access to the catalyst active centers after the deposition of gas phase molecules on it, taking into account the proposed time range of experimental studies. The identification of these stages is directly related to the change in the slope of the curve, which reflects the rate of the reaction. Simultaneously, the reaction rate depends on the amount of free active sites of the catalyst. At the initial stage, the growth rate is approximately constant, which indicates the approximate constancy of the involved active centers. At the second stage of the process, a significant decrease in the growth rate of CNTs is observed, which is associated with the clogging of some of the catalyst active centers with amorphous carbon, which is unable to integrate into the forming nanostructures. During the third stage of the process, a return to the linear growth of CNTs is observed, which indicates a stable growth rate. This occurs because of the establishment of an equilibrium between the processes of formation of amorphous carbon and its removal from the catalyst active sites upon interaction with hydrogen. Notably, the duration of each stage does not have a clear time frame and may vary depending on the synthesis conditions.
The mentioned above key elements of the physicochemical nature of the process under consideration have been proposed to be described using a compact equation: where W0(CH4;H2;T) is function of the slope angle of the CNTs growth first stage on the process control parameters; Kkt(t) is catalyst activity coefficient, accounting the growth rate change over time for stages 2 and 3. The slope of stage 1 of the growth curve strongly depends on the composition of the gas phase and temperature. Studies in the plane of (CH4;H2) parameters were carried out at the same temperature (Fig. 1), and identification of temperature dependencies -for one pair of concentrations of methane and hydrogen. This plan of the analyzed experimental studies allowed to write the function of three variables W0(CH4;H2;T) as the product of the function of two variables WK(CH4;H2) and the function of one variable WT(T): where km is "scaling" factor and is equal to the inverse value of the average initial CNTs growth rate for the "central" experiment, guaranteeing the correct overlap of the dependencies WK(CH4;H2) and WT(T): Expression (2)   The function WK(CH4;H2) was chosen as a two-dimensional cubic polynomial using the co-ordinate descent method (as it was revealed in the course of computational experiments, the quadratic polynomial was not sufficient to describe this dependence): A separate problem in finding the dependence (5) was that the response surface consists of multiple local extrema. Therefore, an additional constraint was imposed on the solution, taking into account the mutual position of the last point on each experimental curve and the hypothetical prolongation of stage 1 to the end of the experiment corresponding to this point in time. Due to the destructive processes occurring at stage 2 when it exists, the hypothetical final point on the calculated curve should not be higher than the real experimental one for all parameter values.
To describe the function Kkt(t), we carried out estimation and analysis of the ratio of the current growth rate We(t) to the initial one for a given set of parameters: where the numerator was calculated on basis of experimental data, and the denominator was estimated on basis of equation (2), taking into account dependencies (3)- (5).
Typical characteristics of the experimental Kkt(t) dependence in the process under study are illustrated by the curve shown in Figure 4. In the first half an hour the peak means some unclear incubation period. However, in some cases it can be absent, and then the curve (6) will be a non-linear decreasing function. We can assume that the peak is always available, but it can be located both in the area of "positive" time (and, accordingly, observed on the graphs) and in the area of "negative" time or it can coincide with the beginning of its counting. Thus, a Cauchy distribution function modified by horizontal and vertical shift parameters was proposed to reproduce the dependence estimated from experiments (6): where γ (scale parameter), t 0 (extremum shift in time parameter), and f 0 (vertical shift parameter) are functions of the process parameters -concentrations of methane, hydrogen, and temperature. Moreover, due to the features of the experimental research plan mentioned above, each of them can be represented according to the principle of equation (2) with its "scaling" coefficients (3).
As computational experiments have shown, the selection of parameters values of equation (7) γ, t 0 , and f 0 for each individual CNTs growth curve provides a large number of alternative solutions that do not differ significantly in accuracy. Therefore, a separate problem was the choice of such solutions set, which in the coordinate plane (CH4;H2) would represent characteristic trends, considering the need for their further approximation, as well as the availability of physical sense for these trends. The found solutions sets f 0 (CH4;H2), γ(CH4;H2) and t 0 (CH4;H2) are presented in Tables 1-3.  The most complicated situation was with the dependence of the vertical shift parameter f 0 (CH4;H2), since its variation had an extremely sensitive effect on the position of the CNTs growth curves. In this case, the effect of the f 0 value on the CNTs growth curve was quite simple: increasing f 0 contributed to more CNTs formation, while decreasing f 0 contributed to less. The f 0 value also affected the curvature of the growth curve: a decrease in f 0 contributed to an increasing of its curvature, i.e. the intensification of the stage 2. In this case, it is possible to single out a conditional diagonal corresponding to the concentration ratio CH4:H2 = 2:1, where the clogging of the catalyst is minimal. That is, we can talk about the optimal ratio of the reactants concentrations as the main criterion for the required dependence f 0 (CH4;H2); deviation from the mentioned proportion 2:1 leads to deterioration of process parameters and complicates its course due to intensification of certain sorption processes.
The solutions set for the scale parameter γ(CH4;H2) ( Table 2), as it was found, can also be connected with the proportion CH4:H2 = 2:1. It is a boundary state, the overcoming of which in the area of increasing methane concentration, as well as decreasing f 0 , contributes to increasing the curvature of the CNTs growth graph and intensification of stage 2. This is associated with high clogging of active centers by amorphous carbon and at the same time weakening the desorption process by hydrogen. However, unlike f 0 , varying the γ value did not change the slope of the growth curve as a whole and did not have such a significant effect on the total number of CNTs.
The solutions set t 0 (CH4;H2) for the horizontal shift parameter (Table 3) identifies the feature of one set of values (CH4;H2), at which the total concentration of active reagents in the gas flow is below a certain critical value, causing the shift of velocity growth peak to a "negative" time, but in fact we observe an extremely inefficient process. In other cases, the growth rate peak can be compared to the time point, which means the absence of the incubation period.
Thus, the sum of (CH4 + H2), which loses its sensitivity limit in the range from 14.3 to 28.6, could be the criterion for the sought dependence t 0 (CH4;H2); for the dependence γ (CH4;H2), criterion was an expression sensitive to the change of sign: = 2 4 − 0.5 (8) and for the dependence of f 0 (CH4;H2), the criterion was an expression of type, which is sensitive only to the deviation from the CH4:H2 = 2:1 proportion: The modulus in expression (9) was added to make it possible not to restrict the potential m values to a set of integers; as subsequently obtained, the optimal m = 0.5.
Approximations of the dependencies t 0 (CH4;H2) and γ(CH4;H2) were sought among the functions of the exponential type, which are sensitive to the sign of the index, and the dependence -in the form of a polynomial. The result of the approximation was the following set of functions: ( 4 ; 2 ) = 25 − 0.1 2.9−5 , (10) 0 ( 4 ; 2 ) = −0.042 3 + 0.1212 2 − 0.272 − 0.0034, (11) 0 ( 4 ; 2 ) = −0.01 11−41.5( 4 + 2 ) . (12) The task of finding and choosing solutions f 0 (T), γ(T) and t 0 (T), in addition to the above mentioned restriction on the physically valid trend availability, should have another restriction related to the matching values possibility with the found dependencies (10)- (12) at 750 ºC. The influence of f 0 , γ, and t 0 values retained the same trends, but given the more pronounced stage 2 at increasing temperature, the calculated CNTs growth curves position was more influenced by parameter γ. At γ values less than 1, the catalyst efficiency and the final CNTs output were significantly reduced, which was just required to describe the process at high temperatures. The found set of solutions is presented in Table 4. Approximation of the dependence f 0 (T) was found in the form of a polynomial, and to approximate the dependence γ(T), the data in Table 4 were prelogarithmic, after which they were also approximated by a polynomial. The data in Table 4 for t 0 corresponded to both the obvious absence of the incubation period and the start of the process at this period end, so it was decided not to find an additional equation for the dependence t 0 (T), but simply to embed the result in function (12). The obtained dependencies f 0 (T) and γ(T) are presented below: 0 ( ) = 0.0002 3 − 0.0004 2 − 0.0137 + 0.0541, (13) ( ) = 0.4572 5 −19.241 4 +322.68 3 −2694.6 2 +11195 −18487 , (14) where t = T/100 -scaling of the temperature axis, made to increase the polynomial coefficients significance at large degrees.
Thus, the proposed mathematical model includes 10 equations containing 38 approximation equation coefficients. It should be noted that the measurements presented in [16] were continuous, therefore, the conditional number of experimental points depended on the experimental curves discretization frequency. The discretization step chosen at the beginning of the work was 0.05 h, which allowed us to interpret the data from [16] as containing 738 conditional experimental points, considering all CNTs growth curves. Figures 5 and 6 show a comparison of our calculations results with the experimental data from [16]. Due to the lack of data in [16] on the experimental studies reproducibility it is difficult to assess the model adequacy, however, the adequacy dispersion was calculated, which was 9.3•10 -4 .  We would like to note that the presented mathematical model is fundamentally different in its structure from the classical kinetic and thermodynamic models, which include directly the equations of physical chemistry and related scientific knowledge areas. The authors of the present work offer the concept of building the mathematical description directly on the experimental data analysis, as well as the various factors influence on the process rate. The approach is based on equations that do not directly describe the mechanism of the process, but take into account the main observed trends. As already mentioned, the main difficulty in finding such equations and estimating their parameters is the solution multiplicity, of which it is necessary to select those that have typical dependence parameter trends and indirectly represent the physicochemical nature of the described phenomena.

Results and discussion
At the same time, approximation of the found sets of solutions can also allow some variability. A clear advantage of the approach is the step-by-step estimation of the numerical equation parameters, which, as a rule, are lacking in classical kinetic and thermodynamic models. For this reason, the need to use optimization algorithms, considering the preliminary analysis of the response surfaces, disappears and is replaced by more efficient approximation methods, such as LSM. This model should not be weaker than classical analogs in terms of prognosis and/or optimizing the process under study. However, a fundamental change in the process conditions, such as the replacement of a catalyst, will require new studies, since there is no universal information in terms of the values of kinetic and thermodynamic parameters.

Conclusion
In this article, we have applied an approach to the development of a mathematical description of chemical processes, which is based on the creation of an approximation equations set, considering a thorough analysis of the features and interactions of the reaction process. Growth of carbon nanotubes was studied in the three-dimensional space of parameters -the methane concentration, the hydrogen concentration in the gas phase, and temperature. The approach assumes a system of equations and dependencies of the rates of various stages of process on its parameters. Based on the results obtained, it can be concluded that this approach can compete with classical methods.