Reaction Molecular Dynamics Study of Combustion Mechanism in Heavy Oil Thermal Recovery

1. Introduction

Karamay heavy oil is not only abundant but also renowned for its unique rheological characteristics. In addition to typical features such as high density, viscosity, and a low freezing point—common in heavy oils from various oil fields—Karamay heavy oil is distinguished by its exceptionally low levels of asphaltene and sulfur. This specific type of heavy oil is highly sought after as a feedstock for producing premium lubricant base oils, resulting in a relatively higher market price compared to other crude oils. Karamay heavy oil reservoirs are typically located at depths of less than 1,000 meters and at lower temperatures, which contribute to their high viscosity under reservoir conditions. Consequently, thermal techniques such as cyclic steam stimulation, steam drive, and in situ combustion are employed to extract this heavy oil. Notably, the viscosity of Karamay heavy oil exhibits less variation with temperature at elevated levels, setting it apart from heavy oils found in other oil fields. While the distinct rheological properties of Karamay heavy oil are believed to be influenced by its chemical composition, further research is necessary to fully understand this phenomenon.

Olalekan [1] and colleagues researched new technologies for heavy oil extraction. Commercial production of heavy oil has heavily relied upon the steam injection method due to its effectiveness in reducing viscosity and improving mobility. However, due to uncertainty of oil price, high energy cost, water requirement, heat loss, and environmental issues, attention is being given to new thermal EOR technologies such as in situ steam generation using thermochemical fluid (TCF) injection. This technology essentially involves downhole generation of heat and pressure from exothermic chemical reactions. respectively. In comparison, from the injection of steam generated at 250°C, the pressure at the inlet of the core was 212 psi and the recovery factor was 71% OOIP. These results therefore confirm the increase in interest in the application of the TCF injection technology in heavy oil production.

Zhao [2] and other studies indicate that as the heating rate increases, the pyrolysis of oil shale occurs at higher temperature ranges. This trend is particularly pronounced with higher oil content. The stability of oil shale during pyrolysis is influenced by both oil content and the pyrolysis atmosphere. Specifically, higher oil content correlates with greater stability in the pyrolysis process. In a nitrogen atmosphere, the pyrolysis interval of oil shale is more concentrated, whereas exposure to air extends the pyrolysis interval and reduces the pyrolysis stability index. Furthermore, an increase in heating rate enhances the release characteristics of the products, which remain largely unaffected by oil content.

These ferroan carbonates are preferentially concentrated at distances ranging from 0.5 to 1.0 meters from sandstone-mudstone contacts. Kinetic modeling effectively reproduces the evolutionary pathways of various diagenetic minerals in the Es4s sandstones, which constitute the upper part of the fourth member of the Eocene Shahejie Formation. More importantly, the modeling results indicate that the dissolution of non-ferroan calcite cement does not improve overall reservoir quality, as it results in the precipitation of more stable ferroan calcite and ankerite cement [3].

Zhang [4] revealed that increased pressure enhances the generation of light components and accelerates the conversion of alkenes to alkanes and aromatics. The release temperature of gases increases under pressure, which promotes the production of alkane gas while simultaneously reducing hydrogen yield. As the pyrolysis temperature rises, the combustion performance of semi-coke deteriorates under low pressure; however, this trend is reversed under high pressure. Finally, a pyrolysis mechanism for the coupling of temperature and pressure in oil shale was proposed.

High temperatures induced various chemical and physical changes that led to the formation of the coking zone. Additionally, the small-aperture channels within the coking zone became concentrated with residual oil as a result of this formation process, which decreased fluidity and made production more challenging [5].

Lab-scale experiments and the analysis of core samples are the primary methods used to investigate compositional changes associated with in-situ combustion [6,7,8,9,10,11].

Furthermore, Figure 1 illustrates that benzaldehyde (C₇H₆O) was detected in the vicinity of the combustion interface (547.8 m). Simultaneously, naphthalene (C₁₀H₈) was identified in samples collected from the unswept area (<547.8 m). Zhao and colleagues [12] employed gas chromatography-mass spectrometry (GC-MS) analysis to confirm the presence of naphthalene in the unswept layers and benzaldehyde at the combustion interface. This finding supports the hypothesis that the pyrolysis and oxidation reactions involve C-O and O-H bonds. Collectively, these results demonstrate the potential for the conversion of benzaldehyde and naphthalene. However, there are limited studies in the existing literature that explore the reaction pathways and networks governing the molecular transition between these two compounds (benzaldehyde and naphthalene).

Enhancing oil recovery can be accomplished by gaining a deeper understanding of the mechanisms involved in the transformation of chemical molecules during the primary reactions associated with heavy oil thermal recovery. Research plays a crucial role in this process. Analyzing complex chemical networks in polymerization, combustion, and chemical engineering has always been essential for elucidating these mechanisms.

Most heavy oil combustion reaction mechanisms operate on a large scale and are primarily macroscopic in nature. The conversion between component molecules mainly involves the number of carbon atoms; however, the reaction pathways and networks between specific molecules have not been thoroughly investigated. By employing molecular dynamics methods to study the reaction processes involved in the transformation of heavy oil combustion, we can gain insights into various reactions occurring during combustion from the perspective of microscopic molecules.

Molecular dynamics methods can simulate networks consisting of thousands of reactions and hundreds of species. However, accurate modeling poses significant experimental and computational challenges. Reactive molecular dynamics (MD) has emerged as an essential tool for understanding these complex networks. In this study, we present a methodology for deriving both a network of reactions and its kinetic parameters from reactive MD simulations.

ReaxFF-MD, a method capable of addressing reactivity in computational calculations, has been proposed as a result of ongoing research developments [13,14]. In carbon/hydrogen/oxygen (C/H/O) fuel systems, such as polystyrene [15], cellulosic materials [16], styrene-butadiene [17], JP-10 [18], polycyclic aromatic hydrocarbons (PAHs) [19], hydrogen [20,21], and various other systems [22,23,24,25,26,27,28], ReaxFF-MD is recognized as a significant supplementary approach to experimental methods.

2. METHODOLOGY

2.1. ReaxFF-MD Reactive Force Field

ReaxFF-MD simulates the formation and breaking of chemical bonds by incorporating the concept of bond order, distinguishing it from non-reactive force fields. Each term, including bond order, is dependent on atomic connectivity. Non-bonded interactions, such as van der Waals and Coulomb forces, are considered between every pair of atoms due to their independent connections. Equation (1) [29] outlines the simulated system utilized by ReaxFF-MD to assess atomic charges through the Electronegativity Equalization Method [30].

E_system=E_bond+E_over+E_angle+E_tors+E_vdWaals+E_coulomb

(1)

In this case, the symbols E_bond, E_over, E_angle, E_tors, E_vdWaals, and E_coulomb stand for bond energy, electrostatic energy, dispersive energy, angle energy, over-coordination energy, and torsional angle energy, respectively. Further evaluation of the force field details is possible [31,32].

2.2. Simulation Details.

To investigate the naphthalene oxidation process, we filled a cubic box with 20 naphthalene (C₁₀H₈) molecules and 240 oxygen (O₂) molecules using Materials Studio 8.0. The volume of the box was calculated based on the initial number of molecules and the target density of 1.00 g/cm³.

We utilized Material Studio 8.0 to populate a cubic box with 20 naphthalene (C₁₀H₈) molecules and 240 oxygen (O₂) molecules to investigate the naphthalene oxidation process. Based on the initial molecular count and the target density of 1.00 g/cm³, we calculated the box volume. The dimensions of the simulation boxes were measured in angstroms (Å). In the subsequent step, we assigned atomic velocities that conformed to a Gaussian distribution at the target temperature. The NVE ensemble, utilizing the Berendsen thermostat, was employed to maintain the system at the desired temperatures for one second. For molecular recognition, we established a C-C bond order cutoff of 0.55 Å, a C-O bond cutoff of 0.65 Å, and a time step of 0.1 fs. To ensure adequate collisions, the simulation temperatures were artificially elevated to 3000 K, 3500 K, and 4000 K. The high temperatures facilitated sufficient collisions between molecules and intermediates, thereby enhancing the simulation process. Within the temperature range of 3000 K to 4000 K, we conducted ten annealing simulations.

Compared to naphthalene, atomic charge analysis indicates that the activating group functions as an electron donor, while the inactivating group acts as an electron acceptor. The effects of substituents on naphthalene are primarily assessed through geometric and vibrational analyses. Furthermore, Natural Bond Orbital (NBO) studies have demonstrated that 2-naphthol esters possess the highest stability energy, which can be attributed to the increased electron density on the oxygen atom. A higher stability energy correlates with increased reactivity of the compound and a smaller band gap [33].

The species was analyzed, and the mechanisms and pathways of the oxidation process were investigated using the LAMMPS package (Large-scale Atomic/Molecular Massively Parallel Simulator) [34,35]. By assessing the number of reaction events that modify the chemical composition and current concentration of the respective reactants, these simulations can be utilized to compute rate constants with software such as the Chemical Trajectory Analyzer (CTY), which was employed in this study to analyze changes in chemical composition within a molecular dynamics (MD) trajectory [36,37].

3. Results and Discussion

A series of NVE-MD simulations of naphthalene oxidation at elevated temperatures (3000 K, 3500 K, and 4000 K) were conducted as described in the section on computational details. Additionally, we employed the CHO-2008 force field to simulate the same system in order to investigate the initiation mechanism and combustion dynamics. During the high-temperature simulations, Figure 2 illustrates the distribution of reactants, products, and intermediate species.

Several reactions, including pyrolysis, occur concurrently with the naphthalene oxidation reaction. This investigation examined both pyrolysis and oxidation reactions. It focused on the effect of temperature on the variation in the quantities of CO₂, CO, H₂O, and other products over time. As the temperature rises to 4000 K, the reaction rate of C₁₀H₈ in O₂ oxidation reaches its peak, and at 60000 fs, the concentration of CO in the reaction exceeds that of CO₂. Additionally, the number of small molecules and the diversity of free radical colors in the oxidation reaction significantly increase at 4000 K.

The three stages of the oxidation process can be identified based on the reaction rate of naphthalene oxidation. In the first stage, naphthalene reacts at a relatively slow rate, ranging from 10,000 fs to 40,000 fs, at temperatures of 3000 K, 3500 K, and 4000 K. Although the production of products is limited at this point, the temperature still influences the slope of the reaction curve. In the second stage, the number of products increases significantly, and the reaction rate of the reactants accelerates. Our simulations revealed the presence of various species, including highly unsaturated compounds, unstable radicals, alkanes, and alkenes. While the type of products remains unchanged in the third stage, the reaction rate does vary. The examination of Figure 2 highlights intriguing topics for further research, such as naphthalene's post-oxidation products and intermediates, their effects on the oxidation reaction, and the transformation of naphthalene into benzaldehyde. At 4000 K, the curve is notably steep due to the abundance of products. We intend to explore the kinetic details of naphthalene's oxidation reaction, particularly at 4000 K, given the active products and reactions involved.

3.1. Detailed Analysis of the Kinetic Mechanisms of the Oxidation Reaction of Naphthalene

The understanding of the complex pyrolysis process relies on the initiation of the MCH oxidation reaction, which determines the subsequent distributions of intermediates and products. Therefore, a comprehensive description of the initial reaction pathway is essential for elucidating the mechanisms of the MCH oxidation reaction.

Figure 2c illustrates the naphthalene oxidation process, which consists of three distinct stages. As previously mentioned, the time axis divisions for these stages remain consistent across various reaction temperature conditions. The rapid oxidation observed in Stage 1 can be attributed to the presence of highly reactive free radicals generated during the reactions [38]. This study identified several active free radicals, including H₂O₂, HO₂, CH₂O, among others. However, we were unable to identify the isomers present in the reaction products, which did not impact the overall objectives of the investigation. By analyzing the data collected at each stage, we were able to determine the reaction rates for the different phases of the process.

Naphthalene undergoes an oxidation reaction that produces carbon dioxide (CO₂), water (H₂O), carbon monoxide (CO), and hydrogen peroxide (HO₂) as the principal products. These four products were analyzed in terms of their reaction types and rates. The highest reaction rate, along with a comparatively high number of reactions, is observed in the reactions listed below.

Naphthalene primarily dissociates through two types of reactions, as illustrated in Table 1: unimolecular decomposition and hydrogen abstraction reactions, both of which involve the dissociation of C-H and C-C bonds. The reactions listed in the table represent those with the highest frequency observed during the simulation process, and they predominantly contribute to the generation of the corresponding products. Naphthalene undergoes an oxidation reaction in which H-O and C-O bonds are formed. During the hydrogen abstraction reaction, the formation of C-O and O-H bonds occurs as part of the main product formation process.

3.2. Intermediate Reactions

Although naphthalene oxidation has been the subject of extensive research, a clear understanding of the intermediate processes has yet to emerge. Comprehending the intermediate mechanisms necessitates a thorough examination of the reactions related to the primary products. The reaction pathways for the main products were studied, and the proportions of reaction times for these products were calculated and documented. The primary forward and reverse reactions in each pathway were discussed, and the key reactions within the overall process were identified.

It is explained how free radicals influence the rate of reaction and synthesis activity. In Figure 2c, carbon dioxide (CO₂) is identified as the most prevalent product resulting from oxidation. The changes in carbon dioxide concentrations throughout the reactions are highlighted. The time-dependent distributions of carbon dioxide are illustrated in Figure 2c, showing a dramatic increase in production between 0 and 40,000 femtoseconds (fs) during the reaction process. During this period, the production of carbon dioxide significantly exceeds its consumption. Figure 3a indicates that reversible reactions are present. Among the various positive reactions, the most frequent is the oxidation reaction of free radicals (CHO₂ + HO → H₂O + CO₂), along with several pyrolysis reactions (C₂O₄ → CO₂ + CO₂). These reactions involve the dissociation of carbon-carbon (C-C) bonds and hydrogen abstraction. After 80,000 fs, the reaction rate of CO₂ approaches equilibrium, at which point the consumption of carbon dioxide surpasses its production. Reversible reactions serve as the primary pathways for the consumption of CO₂.

In the reactions illustrated in the figure, there are 10 positive reactions and 12 negative reactions. Among the positive reactions, the reaction with the highest frequency accounts for 38% (CHO₂ + HO → H₂O + CO₂). In the negative reactions, the reaction with the highest frequency accounts for 34% (H₂O + CO₂ → CHO₂).

In the oxidation reaction of naphthalene, the behavior of carbon monoxide (CO) undergoes significant changes. The production of CO has been continuously increasing during the time interval of 0 to 42,000 femtoseconds (fs), as illustrated in Figure 2c. The reaction with the highest frequency during the CO generation process is the pyrolysis reaction (CO₃ → O₂ + CO), depicted in Figure 3b. Conversely, the production of CO decreases during the period from 42,000 to 89,000 fs, indicating an increase in CO consumption during this stage. The reactions occurring during this interval are predominantly oxidation reactions. A comparison of the curves for CO and carbon dioxide (CO₂) reveals an intersection point, where the production of CO₂ is lower than that of CO within the range of 0 to 62,000 fs. Reason: The revisions enhance clarity, improve technical accuracy, and correct grammatical and punctuation errors while maintaining the original meaning of the text.

There are 16 positive reactions and 5 reverse reactions. Among the positive reactions, the one with the highest frequency accounts for 47% (CO₃ → O₂ + CO). In the case of reverse reactions, the reaction with the highest frequency constitutes 69% (O₂ + CO → CO₃).

In Figure 2c, there is a significant increase in H₂O production during the initial stage of 0-30,000 fs. In the subsequent stages, H₂O production remains relatively stable. In the positive reaction channel depicted in Figure 3c, the proportion of reaction counts reached 33% (H₂O + HO + HO → H₂O + H₂O + O). Conversely, the reverse reaction accounts for 42% of the negative reaction channels. Overall, the ratio of positive reactions producing H₂O is higher than that of the negative reaction channels. This disparity is one of the reasons for the relative stability in H₂O production. Additionally, HO radicals participate in numerous reactions.

The generation curve data for the HO₂ reaction is illustrated in the locally enlarged image in Figure 2c. After 35,000 fs, there is a significant increase in HO₂ production. An analysis of the reaction channels presented in Figure 3d reveals that HO₂ is involved in numerous oxidation reactions, as well as in the thermal decomposition of hydrocarbons and the dissociation of certain free radicals. These findings indicate that HO₂ exhibits strong reactivity and can facilitate various reactions. Within the reaction channel for HO₂, there are 11 forward reactions and 10 reverse reactions. The most prevalent reaction is the oxidation of hydrogen peroxide by oxygen (O₂ + H₂O₂ → HO₂ + HO₂), with the reverse reaction occurring at a frequency that accounts for 27% of the total.

3.3. Detailed Reaction Map for Naphthalene Oxidation

Figure 4 illustrates reaction networks that have been suggested for the oxidation of naphthalene(C₁₀H₈) based on the outcomes of our simulation.

The impact of isomers was overlooked during the research on reaction networks. The reaction rate is high, and the reaction time is short, based on the simulation conditions, which include elevated temperatures. Within the reaction network, isomers are believed to have minimal influence. These networks provide a comprehensive explanation of the process that converts a C₁₀H₈ molecule into C₇H₆O. Four potential pathways are proposed for the formation of C₇H₆O.

The P₁ pathway is the most straightforward of the four reaction pathways, as illustrated in Figure 4. The compound C₇H₄O is essential to the P₁ reaction pathway and is involved in the P_1-2, P_1-4, and P_1-5 reactions, among others. Additionally, C₇H₅O₂ is critical to the P_1-3, P_1-6, and P_2-5 reaction pathways, which presents a similar challenge. C₇H₅O₂ is produced at the intersection of the P₁ and P₂ reaction pathways.

The P₂ reaction pathway incorporates the products from both the P₃ and P₁ reaction pathways. C₁₀H₇O₃ serves as the initial reactant in the P₂ reaction pathway and plays a crucial role in the subsequent reactions. The P₃ reaction pathway is relatively complex, encompassing numerous products, including eight smaller reaction pathways. The product that lies at the intersection of the P₃ and P₂ reaction pathways is C₁₀H₇O₃.

The P₄ reaction pathway is relatively short and is the shortest among all reaction pathways. It exhibits the highest conversion efficiency and a greater reaction rate. The P₁, P₂, and P₃ reaction pathways encompass various types of reactions, including oxidation reactions, pyrolysis reactions, and hydrogen transfer reactions. In contrast, the P₄ reaction pathway consists solely of pyrolysis reactions. The specific reactions and their corresponding rates involved in these four reaction pathways will be detailed in Table 2, Table 3 and Table 4.

Table 2 through Table 4 present the intermediate reactions for each reaction pathway, along with their corresponding reaction rate constants. The reaction map indicates that a total of 20 intermediate reactions are involved. There are three categories of statistical reactions: hydrogen abstraction, pyrolysis, and oxidation. In some of these reactions, hydroxyl (OH) and hydrogen (H) radicals play a significant role.

Based on the reaction network map and Table 2, we will discuss the characteristics of the P1 reaction pathway. This pathway consists of six reactions. The initial oxidation reaction, represented by the equation O₂ + C₁₀H₈ → C₁₀H₈O₂, exhibits the lowest reaction rate. The intermediate product, C₇H₄O, is formed through two reactions, P_1-2 and P_1-5, both of which involve the pyrolysis of C-C bond dissociation. Additionally, two reactions contribute to the formation of the product C₇H₅O₂, specifically P_1-3 and P_2-5. The hydrogen elimination reaction occurs in the P_1-4 pathway, resulting in the production of C₉H₅O₂, as observed.

As a node in the reaction network, C₇H₅O₂ has an alternative generation pathway, as illustrated in Figure 4. According to Table 3, three of the reactions are classified as pyrolysis reactions. Additionally, there is one oxidation reaction and one activation polymerization reaction. The oxidation reaction exhibits the lowest rate among these processes. The end products of the pyrolysis reaction (C₈H₇O₂ → C₇H₅ + CH₂O₂), which results from the dissociation of C-C bonds, are C₇H₅ and CH₂O₂. This particular reaction has the highest reaction rate. A relatively high reaction rate is also observed when the free radical OH is involved.

Compared to the previously mentioned reaction pathway, the P₃ reaction pathway is more complex, involving a total of nine reactions. Among these, the proportion of pyrolysis reactions is relatively high. The oxidation reaction exhibits the slowest rate among all the reactions. When examining the hydrogenation (C₇H₅O₃ + H → C₇H₆O₃) and dehydrogenation (C₁₀H₈O₃ → H + C₁₀H₇O₃) reactions, it is evident that the rate of the hydrogenation reaction is significantly higher than that of the dehydrogenation reaction. Conversely, the rate of the deoxidation reaction (C₇H₆O₃ → C₇H₆O + O₂) is greater than that of the oxidation reaction (O₂ + C₁₀H₈O → C₁₀H₈O₃).

As illustrated in Figure 4, the transition from C₁₀H₈ to C₇H₆O involves two critical products: C₇H₅O₂ and C₇H₅O₃. A comparison of the two reactions involving C₇H₅O₂ reveals that the oxidation reaction (C₇H₅ + O₂ → C₇H₅O₂) has a relatively low rate, as shown in Table 4. In contrast, the formation of product C₇H₅O₃ occurs at a higher rate during the pyrolysis reaction, which involves the participation of free radical OH. Among the two products mentioned, C₇H₅O₂ reacts with free radical OH to produce C₇H₆O₃, while C₇H₅O₃ is formed through a hydrogenation reaction. Ultimately, C₇H₆O₃ yields the target product, C₇H₆O, via a pyrolysis reaction (C₇H₆O₃ → C₇H₆O + O₂).

4. Conclusions

The ReaxFF molecular dynamics (MD) method was employed to investigate the distribution of products during the oxidative pyrolysis of naphthalene. Additionally, the reaction pathway and characteristics of oxidative pyrolysis were analyzed using the CTY method. The reaction products and rates exhibit distinct characteristics at varying temperatures. The reaction pathway network for the conversion of C₁₀H₈ to C₇H₆O primarily comprises three pathways, which include ring opening, carbon-carbon (C-C) bond cleavage, hydrogen transfer, and other related processes.

• Naphthalene undergoes oxidation at elevated temperatures using NVE-MD (microcanonical ensemble molecular dynamics) with a reactive force field. To obtain a comprehensive understanding of the oxidation process, simulations are conducted on a model system consisting of 20 naphthalene molecules. The primary species identified in the current simulations include CO₂, CO, H₂O, HO₂, and others.
• The reaction rates of products in oxidation and pyrolysis reactions vary at different time stages. After the intersection of CO₂ and CO, the yield of CO₂ surpasses that of CO. This phenomenon can be attributed to insufficient collision contact in the early stages and more favorable reaction conditions in the later stages. In the reaction network from C₁₀H₈ to C₇H₆O, multiple types of reactions were identified across four reaction pathways, including dehydrogenation, hydrogenation, and small molecule activation reactions.
• The reaction network from C₁₀H₈ to C₇H₆O was examined. The findings of this investigation will enhance our understanding of the process. Among the reactions in the network, the hydrogenation reaction exhibits the highest reaction rate. Identifying the reaction network through conventional experiments is challenging due to the complexity of the reactions involved. In the study of oxidative pyrolysis reactions, computational simulations offer a novel research approach.