1. Introduction
Fischer–Tropsch synthesis (FTS) is a process used to produce a vast range of organic compounds. Discovered nearly a century ago, FTS reacts H
2 and CO over a metal catalyst to produce C
1-C
70-chain compounds, including gasoline and alternate materials [
1,
2]. Thus, it is often used to generate fuels, monomers for common polymers, rubber, etc. Despite the long history of research and industrial use of FTS, it is still unknown how this process works on a microscopic level [
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
16,
17,
18,
19,
20,
21,
22]. There are at least three proposed mechanisms, namely, the carbide originally proposed by Fischer and Tropsch, CO insertion, and formate mechanisms. Yet, none of these can fully describe the entire FTS process, mainly because of the diversity of materials generated by different active metals [
1].
To further comprehend FTS, many have attempted to use mathematical modeling as a useful means of describing the process [
23,
24,
25,
26,
27,
28]. Density functional theory (DFT) calculations, paired with experimental procedures, have provided further insight into the FTS mechanism [
18,
19,
20,
21]. These models generally fall into one of two categories—probabilistic models and kinetic models. Probabilistic models attempt to describe the propagation of the carbon chain to generate the products observed [
23,
24,
25]. Kinetic models focus more specifically on the rate of the reaction, determining a mechanism to describe the reaction as a whole [
26,
27,
28]. While these models have significant application for understanding FTS, they do not exist without limitation. Kinetic models depend on the mechanism of the reaction, which, as stated earlier, is not fully understood due to the diverse array of products. Furthermore, they become quite mathematically complex. Probabilistic models focus on the likelihood of a specific carbon chain propagating or terminating. This type of modeling is limited by the monomer used for propagation, which is not fully understood.
This research sought to create a model that describes the rate of change of hydrogen and carbon monoxide as the reaction progresses. Using a series of industrially relevant cobalt-catalyst runs with varying supports, a system of differential equations was created to model the unreacted syngas in terms of the products generated. Each chemical product output term (such as moles of methane, water, carbon dioxide, etc.) was scaled by an impact parameter, a value in units of hr-1. The model was then tested using sensitivity analysis (SA) to determine which impact parameter and, consequently, which associated chemical product term altered the solution of the model the most. From the results of this analysis, conclusions were drawn about the chemical implications, which may then be used to further examine catalytic runs.
2. Computational Modeling
The computational model was established using chemical engineering process principles from the stirred-tank slurry reactor and basic chemical kinetics. The rate of H
2 and CO coming out of the reactor (unreacted syngas) was modeled using ordinary first-order differential equations (ODEs) [
29]. More specifically, the rate of change of H
2 in the reaction can be described by the hydrogen flow into the reactor minus everything coming out of the reactor (H
2, water, methane, hydrocarbon products). All the products (excluding unreacted hydrogen) were scaled by an “impact parameter”. This model was applied to CO as well to produce two ODEs as follows:
where Q
H2 and Q
CO are input rates of H
2 and CO in moles per hour, respectively; q
H2 and q
CO are the output rates for H
2 and CO in units of hr-1, respectively; MH
2(t), MCO(t), MH
2O(t), MCO
2(t), MCH
4(t), MC
2C
4(t), and Mliq(t) are the time-dependent functions for the moles of H
2, CO, CO
2, H
2O, CH
4, C
2–C
4 products, and C
5+ (or liquid) products, respectively; and Z1–10 are the impact parameters in units of hr-1. It is important to note that MH
2(t) and MCO(t) are the unknown functions for this system. All other mole functions were interpolated from the experimental data (see
Supplementary Materials for data), but these two are the functions for which the system is solved. The MH
2(t) and MCO(t) solutions of the ODE system were compared to the experimentally observed values of these functions (see
Supplementary Materials) and the resulting sum of squared error was the “response” for the sensitivity analysis. The * symbol at the end of the hydrogen equation is to note that the bracketed term was added later, thus explaining the variance in sequence of impact parameters. Originally, the carbon dioxide term was not included in the hydrogen equation but was later added to explore the possibility of a water–gas shift (WGS).
Due to the innovative nature of the differential equations modeling the reactor, the solver used to evaluate the solutions to this model was the Livermore Solver for Ordinary Differential Equations (LSODA). LSODA performs a test for numerical stiffness as the time steps evolve. A stiff differential equation is one which exhibits solutions with different time scales. For example, stiffness may occur in chemical modeling when simulating reactions with very different reaction rate constants. To solve stiff equations numerically, one must use an implicit solver that requires time-consuming nonlinear equation solutions (e.g., Newton’s method for nonlinear systems) for each time step. The LSODA solver primarily uses a non-stiff solver for speed but dynamically checks the solution to determine if a stiff method is required and applies an implicit solver for those time steps. Since the values of the impact parameters were unknown, this approach provided the most flexibility in the numerical solutions. Initial testing of the solution showed good qualitative agreement with the data, indicating that the LSODA solver worked properly with this model. Once the solution was established, sensitivity analysis (SA) was performed. SA is a branch of computational methods that provide insight into how much the total change in the output of a model can be divided up and allocated to the “input” variables/parameters of the model [
30]. In other words, SA measures the effect on the output of the model that results from changing the input factors. While there are two types of SA, local and global, local was not used in this work since it involves only measures the change in output at a specified point. Global SA measures the output across a specified range for the input parameters, which was more practical in analyzing this model [
30].
While there are many methods for global SA, the Sobol method and Morris method [
31,
32,
33] were of primary interest for application to the model. These are commonly used for such systems and they are both well validated in multiple applications [
32,
34,
35]. The Sobol method is a variance-based global SA method that generates sensitivity indices from ratios of variance terms [
33]. The limitation of Sobol is that calculating the variances is very complex mathematically, and the method requires a large number of samples (often more than 10,000) to counteract any error from negative indices. Furthermore, the primary function of the Sobol method is to analyze the interactions by looking at the second- and higher-order indices, but since the primary focus of this work was the impact of a factor on the model output, the Morris method was used as the SA tool for this preliminary study.
The Morris method is known as a one-at-a-time (OAT) global SA method [
31]. This method measures the impact of each input factor on the model by changing only one input factor at a time and observing the change in the model output. Specifically, Morris generates a set of elementary effects (EEs) for each input factor across the given range [
31]. An EE is generated by starting with the model evaluated at a specific set of input factors. Then, the model is evaluated again with only one factor increased or decreased by a specified amount, D. The unchanged model output is subtracted from the altered model output, the absolute value is taken for the difference, and this absolute value is divided by D, as follows:
In the equation above, Y represents the model, X
1…i is the set of parameters, and i = {1, 2, …, n}, n being the number of parameters in the model. A set of these EEs using multiple values for D is generated for each input parameter across the given range. Since each set, known as the distribution of EEs, can be quite large, they are sampled. The original Morris sampling technique was a type of random sampling that did not use a high number of samples (10–50), but Campolongo et al. discovered a better sampling technique [
31]. In this new technique, 500–1000 original samples are taken from the distribution of EEs. These are then analyzed to determine the optimal 10–50 samples (those that most closely describe the set as a whole). The optimal samples, or optimal trajectories, create a new set that is much smaller and easier to analyze statistically. The average, µ*, and the standard deviation, σ, are taken for each optimal sample set. µ* represents the impact of the specific parameter on the model, and σ represents the interaction between a parameter and the other input factors [
31]. The greater the magnitude for these values, the greater the impact (µ*) and interaction (σ) for the corresponding input parameter. Since Morris is a relatively simple method that focuses more on the impact and less on the interactions between parameters, it was chosen as the global SA method for this work. It is important to note that, when performing SA on this model, the impact parameters were the input parameters. While each impact parameter has a corresponding output product (e.g., k
1 corresponds to the hydrogen term for H
2O output), there is not a necessary relationship between the degree of impact and the product itself. However, this work hypothesizes that there is at least a small direct correlation between the magnitude of impact for an impact parameter and the amount of corresponding product produced.
3. Results and Discussion
To ensure the effectiveness of this modeling, the mass closure for these runs had to be accounted for. Since this work was to understand, specifically, how the reactants affected each product through the sensitivity analysis, our entire process would collapse and be non-functional if mass closure was not accounted for. To do this, the mass flow controllers were calibrated beforehand, the flow in and out was measured regularly, and all liquids pulled from the system were weighed and accounted for. The difference between output and input displayed almost 96% closure for mass balance for all our work.
Using the data collected from the experimental runs, the mole functions were interpolated, the differential equations were solved, and the solution was then tested using Morris method SA. Initially, the DEs were loosely coupled, meaning that the solutions for MH
2(t) and MCO(t) could be somewhat dependent on one another and the input parameters (Z1–Z9) could potentially interact. Each set of experimental data was analyzed nine times with Morris, on varying ranges from 0.1–50 to 0.1–500. The number of samples generated varied between 500 and 1000, and the optimal trajectories varied between 0 and 50, as described by Campolongo et al. [
31]. The µ* and σ values were collected from the Morris output and the average for all nine analyses was taken. The µ* and σ averages for each data set were compared to the other data sets with the same catalyst–support combination and then compared to the values for the other category of data to discern any variance between the two types of supports.
Next, the two DEs were uncoupled. The CO
2 term for the hydrogen equation, along with k
10, was added to explore the possibility of WGS occurring in the 0.5 wt% Pt, 25 wt% Co/Al
2O
3 catalyst. Both equations were analyzed using Morris, following the same procedure as described above. The µ* and σ values were recorded for both the H
2 and the CO equations. These were then compared, primarily the µ* values for k
5 and k
10, the CO
2 impact parameters for the CO and H
2, respectively. In theory, should µ* for k
10 consistently be greater than that of k5, it could be indicative of hydrogen having a more substantial role in producing CO
2. This could only happen through WGS, and thus if Z10 > Z5, it could indicate the possibility of WGS. Otherwise, if Z5 > Z10, this would most likely perpetuate the prevailing theory, primarily for cobalt catalyst FTS, that the Boudouard reaction is producing the CO
2 [
36,
37].
Listed below are the data collected from the outlined modeling and sensitivity analysis procedure. As a key,
Table 1 is the list of the five product categories and the corresponding impact parameters for the two ODEs in the model. For example, Z
4 is the impact parameter for water in the carbon monoxide equation.
Recall that the µ* value for each parameter indicates how much that parameter influences the model. In other words, changing the impact parameter with the highest µ* value will alter the output of the model more than any of the other parameters. The µ* value does not indicate any level of interaction among the parameters, only how significantly that parameter affects the model. To examine the interactions between the parameters effectively, Sobol method SA should be used.
Table 2 and
Table 3 show the data collected from Morris method SA performed on the coupled DEs.
From these data, the order of impact for the nine impact parameters can be observed for both supports. For 20 wt% Co/SiO
2, the order of impact is Z
7 > Z
3 > Z
4 > Z
1 > Z
6 > Z
2 > Z
8 > Z
9 > Z
5. In contrast, the order for 0.5 wt% Pt, 25 wt% Co/Al
2O
3 is Z
7 > Z
3 > Z
4 > Z
1 > Z
6 > Z
2 > Z
8 > Z
5 > Z
9. Although a subtle difference, the order of k
9 and k
5 is inverted, indicating that k
5 had a greater impact on the model with the cobalt–alumina catalyst than the cobalt–silica. These orders of impact held constant across all three sets of data for each type of support, as can be seen in
Table 2 and
Table 3.
Table 4 lists the averages for each category side-by-side in order to compare the two supports more closely. From the hypothesis, this variation in the order of the parameter impact indicates that CO
2, the corresponding product for Z
5, is more impactful on the alumina support than the silica.
Furthermore, based on the proposed correlation between the products and corresponding parameters, alumina produces a greater amount of carbon dioxide than does silica. To verify this hypothesis, the original data from the data sets were examined, and the theory was confirmed. As
Table 5 shows, the cobalt–alumina catalysts did in fact produce consistently more carbon dioxide than the cobalt–silica.
Due to this CO
2 production difference, the two DEs were separated and examined individually using the Morris method. Only the Category 2 data were explored since the CO
2 parameter was more impactful for the cobalt–alumina support. The additional carbon dioxide must come from an additional source. According to the hypothesis, that source could be either water–gas shift (WGS) or the Boudouard reaction. In order to account for any CO
2 produced as a result of hydrogen interaction, thus allowing for the possibility of WGS, the parameter Z
10 was added to the hydrogen equation along with the CO
2 mole function.
Table 6 and
Table 7 show the data collected from Morris method SA for the separated equations.
The parameters Z5 and Z10 were of primary interest. As is evident from the data, the impact of k5 was invariably greater in magnitude for the CO equation than the impact of Z10 for the H2 equation across all three data sets in Category 2. Since Z5 > Z10, the likelihood of carbon dioxide being produced from WGS is very minimal. The results indicate that the extra CO2 is a product of the Boudouard reaction.
5. Conclusions
In conclusion, the mathematical model created in this work accurately described the FTS process in terms of rate of change of syngas. When the model was analyzed using Morris method SA, a significant difference was observed between cobalt–silica and cobalt–platinum–alumina catalysts, namely a difference in the effect of the impact parameters k
5 and k
9. Based on the initial hypothesis, this difference indicated that the cobalt–platinum–alumina generated more carbon dioxide than did the cobalt–silica. This was verified by the raw data as well. These conclusions are highly specific to the two catalyst–support types in question, and thus more general conclusions must come from further research. Thus, the equations were separated, and the cobalt–platinum–alumina catalyst data were reexamined, now with an additional CO
2 term on the hydrogen equation to test for the possibility of WGS. The observed results showed that the carbon monoxide contribution to CO
2 was consistently higher than the hydrogen contribution. Since WGS is the only pathway in which hydrogen produces CO
2 during FTS, WGS was excluded as a source of CO
2 production. Hence, the Boudouard reaction was determined to be the more probable source of additional carbon dioxide. Though the Boudouard reaction yields a high amount of carbon on the surface, to the best of our knowledge, there is no direct evidence that this reaction leads directly and solely to inert graphitic carbon (coke). Moreover, one of the main mechanistic routes proposed for FT is the carbide mechanism, where CO completely dissociates before reacting with hydrogen. In order to understand this rate of formation for C, evidence of more than one type of carbon on the catalyst surface, an active phase and an unreactive phase, has been observed on multiple catalysts such as cobalt [
40] and Fe [
41].
Further indications using sensitivity analysis show that both the H
2 and CO equations displaying Z
3 and Z
7 obtained the highest value, indicating the liquids are most sensitive to any deviations brought on by H
2 and CO. Though expected, these details could further indicate the importance for a specific balance in the FTS system between the two reactants, suggesting the difficulty in tuning a specific catalyst for C
5+ products as a whole. Additionally, the C
5+ products are more sensitive in the CO reaction than the H
2, which implies that the balance for the FTS system lies more heavily in the CO than for H
2, primarily describing the importance of CO, possibly due to the metal’s ability to back donation [
42] in the FTS process.
While this research is preliminary, the model and analysis produced some excellent descriptive results for FTS. Although this work only examined cobalt catalyst data, future work could include the incorporation of multiple catalysts and supports. Impact parameters could potentially be optimized in future work to find the ideal values for a given catalyst. The model can be adapted to include terms for olefins and paraffins, and the liquid (C5+) term could be broken down into individual terms for each number of carbons per chain (e.g., C5, C6, C7, etc.). The equations could be changed to predict specific product amounts based on initial conditions such as catalyst, support, temperature, pressure, flowrate, etc. Ultimately, this type of FTS modeling not only describes the process in various and versatile ways but can also be adapted and altered to potentially produce a predictive model.