A method and system for construction verification of a vacuum residue model
Patent Information
- Application Number
- CN202611023108.1
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-07-10
- Publication Date
- 2026-10-20
AI Technical Summary
然而,现有研究多集中于煤的热解或单一的沥青模型,针对减压渣油这一复杂体系的平均分子模型构建及其热解机理的系统性研究仍十分匮乏
本申请通过系统分析减压渣油的组成与性质,构建了精确的C57H74OS平均分子模型,填补了减压渣油领域缺乏精确分子模型的技术空白。还提供了从宏观密度到微观键级/电荷的多尺度验证方法,确保了模型的准确性。同时,基于验证后的模型进行ReaxFF-MD模拟,阐明了乙烯和氢气的生成路径差异,揭示了杂原子(O、S)从稳定杂环转化为H2O、H2S等小分子的机理,为减压渣油的资源化利用和工艺优化提供了理论依据。因此,本申请的建模与验证方法不仅适用于减压渣油,也为其他复杂油品的分子模拟研究提供了可借鉴的技术路径。
Smart Images

Figure FT_1 
Figure FT_2 
Figure FT_3
Abstract
Description
Technical Field
[0001] This application relates to the field of residual oil modeling, and in particular to a method and system for constructing and verifying depressurized residual oil models. Background Technology
[0002] Vacuum residue is a heavy residue after crude oil distillation, with an extremely complex composition. It is currently commonly classified into four components: saturated components, aromatic components, gums, and asphaltenes. The aromaticity, from lowest to highest, is: saturated components, aromatic components, gums, and asphaltenes. While this four-component classification effectively characterizes the macroscopic composition of the residue, the resulting molecular structure is complex, and the proportions of each component are difficult to determine, leading to significant modeling challenges and making it difficult to accurately analyze its reaction behavior at the molecular level.
[0003] Currently, the combined use of density functional theory quantum computing and ReaxFF reactive force field molecular dynamics has become the mainstream approach for studying the pyrolysis mechanism of vacuum residue, enabling systematic analysis of chemical bond breaking, free radical recombination, and intermediate evolution during the pyrolysis process. However, existing research largely focuses on coal pyrolysis or single bitumen models, and systematic studies on the construction of average molecular models and the pyrolysis mechanism of the complex vacuum residue system remain scarce. Summary of the Invention
[0004] This application provides a method and system for constructing and verifying a vacuum residue model, thereby filling the gap in the prior art where there is a lack of accurate molecular models for vacuum residue.
[0005] The first aspect of this application provides a method for constructing and validating a vacuum residue model, comprising the following steps: Step 1: Obtain experimental data of vacuum residue, including elemental analysis data, number-average molecular weight data, and infrared spectral data; Step 2: Construct an average molecular structure model of vacuum residue based on the experimental data; Step 3: Perform density verification of the average molecular structure model using molecular dynamics simulations, including performing energy minimization, annealing, and pressure cycling on the system in sequence, and calculating its density under isothermal and isobaric ensembles, and comparing the calculated results with the experimental density for verification. Step 4: Perform quantum chemical calculations to verify the constructed and validated vacuum residue model, including calculating its infrared spectrum, Mayer bond order and CM5 charge, and compare and verify the results with experimental infrared spectra and reaction force field simulations. Step 5: Based on the validated vacuum residue model, construct a pyrolysis reaction model and use reaction force field molecular dynamics simulation to study its pyrolysis behavior.
[0006] Optionally, the construction of the average molecular structure model of vacuum residue in step 2 specifically includes: determining the molecular formula of vacuum residue as C based on elemental analysis data and number-average molecular weight data. 57 H 74 OS; Based on unsaturation calculations and infrared spectral analysis, the molecular structure is determined to be based on a fused aromatic ring, a cyclohexane ring, and oxygen-containing and sulfur-containing heterocycles as the core framework, with alkyl side chains attached; Based on the assumption that each cycloalkane ring and benzene ring is connected to an alkyl side chain and there are no multiple branches, two-dimensional and three-dimensional average molecular structure models are constructed.
[0007] Optionally, the density verification in step 3 specifically involves: minimizing the energy of the system using the conjugate gradient method, then sequentially performing heating-cooling annealing and pressure cycling operations under an isothermal and isobaric ensemble, and finally performing isothermal and isobaric ensemble simulations at 300K and 1 bar to calculate the system density and compare it with the experimental density. If the relative error is within the acceptable range of molecular dynamics simulation, the model passes the density verification.
[0008] Optionally, the quantum chemical calculation verification in step 4 specifically includes: calculating the vibrational frequency of the constructed model using density functional theory, correcting the calculated wavenumber, comparing the peak position and peak intensity of the corrected calculated infrared spectrum with the experimental infrared spectrum; calculating the Mayer bond order and CM5 charge of the model, comparing them with the bond order and charge obtained from the reaction force field simulation calculation, and calculating the mean absolute error and coefficient of determination to verify the reliability of the force field parameters.
[0009] Optionally, the correction of the calculated wavenumber specifically includes: using a correction factor of 0.9525 for antisymmetric stretching vibration in the high-frequency region, using a correction factor of 0.9483 for symmetric stretching vibration in the high-frequency region, and using a correction factor of 0.9793 for vibration in the low-frequency region.
[0010] Optionally, step 5 specifically includes: constructing a periodic pyrolysis reaction model containing multiple molecules based on the validated vacuum residue model; performing molecular dynamics simulations of the reaction force field at different heating rates to analyze the residual mass of reactants, the product formation rules, and the evolution trends of key products ethylene and hydrogen during pyrolysis; and dividing the pyrolysis process into three stages—inert heating, rapid cracking, and slow reaction—based on the simulation results, and analyzing the influence of different heating rates on product distribution.
[0011] Optionally, the different heating rates include 10K / ps, 12K / ps and 16K / ps, with a final heating temperature of 3500K. The optimal simulated temperature range for vacuum residue pyrolysis is determined to be delayed as the heating rate increases, with ethylene yield remaining at a relatively high level in the range of 2900K to 3300K.
[0012] Optionally, the process also includes step 6: based on the results of molecular dynamics simulation of the pyrolysis reaction, analyzing the structural characteristics of vacuum residue, wherein the structural characteristics include: high carbon content, abundant alkyl side chains, high proportion of aliphatic structure residue, heteroatoms existing in the form of aromatic heterocycles in furan and thiophene rings, and the molecular skeleton being mainly composed of polycyclic aromatic hydrocarbons and cyclohexane structures.
[0013] A second aspect of this application also provides a system for constructing and validating a vacuum residue model, comprising: The data acquisition module acquires experimental data of vacuum residue, including elemental analysis data, number-average molecular weight data, and infrared spectral data. The model building module is used to construct an average molecular structure model of vacuum residue based on the experimental data. The density verification module is used to verify the density of the average molecular structure model through molecular dynamics simulation. The quantum chemical verification module is used to verify the average molecular structure model by infrared spectroscopy, Mayer bond order and CM5 charge through quantum chemical calculations. The pyrolysis simulation module is used to build pyrolysis reaction models based on validated models and perform molecular dynamics simulations of reaction force fields.
[0014] Beneficial effects: This application constructs a precise C-type structure by systematically analyzing the composition and properties of vacuum residue. 57 H 74 The OS average molecular model fills the technological gap in the field of vacuum residue, which lacks accurate molecular models. It also provides a multi-scale verification method from macroscopic density to microscopic bond / charge levels, ensuring the model's accuracy. Furthermore, ReaxFF-MD simulations based on the verified model elucidate the differences in the formation pathways of ethylene and hydrogen, revealing the mechanism by which heteroatoms (O, S) transform from stable heterocycles into small molecules such as H2O and H2S, providing a theoretical basis for the resource utilization and process optimization of vacuum residue. Therefore, the modeling and verification methods in this application are not only applicable to vacuum residue but also provide a referable technical path for molecular simulation studies of other complex oil products. Attached Figure Description
[0015] Figure 1 The infrared spectrum of vacuum residue in one embodiment is shown. Figure 2 This is a planar single-molecule structure model of a vacuum residue sample in one embodiment; Figure 3 This is a stereoscopic single-molecule structure model of a vacuum residue sample in one embodiment; Figure 4 This is a schematic diagram of the density distribution of vacuum residue in the last 50 ps of one embodiment; Figure 5 C5H in one embodiment 12 C 10 H 22 C 15 H 32 C 20 H 42 Calculated infrared spectrum; Figure 6 This is a magnified view of the high wavenumber region in one embodiment; Figure 7 This is a magnified view of the low wavenumber region in one embodiment; Figure 8 Here are the experimental and calculated infrared spectra from one embodiment; Figure 9 This is a schematic diagram of key-level comparison in one embodiment; Figure 10 This is a schematic diagram of charge comparison in one embodiment; Figure 11 This is a schematic diagram of the proportions of each component at 10K / ps in one embodiment; Figure 12 This is a schematic diagram of the proportion of each component at 12K / ps in one embodiment; Figure 13 This is a schematic diagram of the proportion of each component at 16K / ps in one embodiment; Figure 14 This is a thermogravimetric analysis comparing simulation with actual results in one embodiment; Figure 15 This is a schematic diagram showing the amount of C2H4 at different heating rates in one embodiment. Figure 16 This is a schematic diagram showing the amount of H2 at different heating rates in one embodiment. Figure 17 This is a schematic diagram of the product path analysis of vacuum residue in one embodiment; Figure 18 This is a schematic diagram illustrating the conditions of vacuum residue oil at different simulated temperatures in one embodiment. Detailed Implementation
[0016] To make the objectives, technical solutions, and advantages of this application clearer, the following detailed description is provided in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative and not intended to limit the scope of this application.
[0017] This application provides a method for constructing and validating a vacuum residue oil model, comprising the following steps: I. Construction of the vacuum residue model: Step S1: Determine the average molecular formula. The mass fractions of C, H, O, N, and S in the vacuum residue were determined using an elemental analyzer; the number-average molecular weight M of the vacuum residue was determined using gel permeation chromatography. n The number of atoms in each atom was calculated based on the mass fraction and number-average molecular weight, and then corrected by combining the functional group types determined by Fourier transform infrared spectroscopy to obtain the average molecular formula of vacuum residue.
[0018] Step S2: Construct an average molecular structure model of vacuum residue. Based on unsaturation calculations, infrared spectral analysis results, and core principles of molecular structure, determine the ring structure combination, side chain length, and position of the molecule.
[0019] Step S3: Generate a 3D structure file. A 3D model of the molecule is constructed using molecular simulation software, exported as an input file for quantum chemistry calculation software (such as ORCA), and then structural optimization is performed.
[0020] II. Validation of the vacuum residue model: Step S4: Density Verification. The constructed vacuum residue model is placed in a simulation chamber, and an NPT ensemble simulation is performed in Lammps molecular dynamics software using the OPLS-AA force field. The density of the simulated system is calculated and compared with the experimentally measured density. If the error is within an acceptable range, the model is considered reasonable.
[0021] Step S5: Infrared Spectroscopy Verification. The constructed vacuum residue model undergoes DFT frequency calculations. After correction by a correction factor, the calculated infrared spectrum is obtained and compared with the experimental infrared spectrum. The accuracy of the model is further verified by the matching degree of peak position and peak shape.
[0022] III. Pyrolysis Reaction Simulation: Step S6: Construction and Validation of the Pyrolysis Reaction Model. The validated model was used for ReaxFF-MD reaction molecular dynamics simulation. The Mayer bond order calculated by the ReaxFF force field was compared with that calculated by DFT, and the CM5 charge calculated by ReaxFF-QEQ was compared with that calculated by DFT to verify the reliability of the force field parameters.
[0023] Step S7: Pyrolysis Mechanism Analysis. ReaxFF-MD pyrolysis simulations were performed on the vacuum residue model at different heating rates to analyze the formation patterns and reaction pathways of key products (ethylene and hydrogen).
[0024] Each step is explained in detail below: I. Construction of the vacuum residue model: Step S1: Determine the average molecular formula: Number-average molecular weight determination: The average molecular formula was obtained by elemental analysis and gel permeation chromatography (GPC). The GPC used was an Agilent GPC 50 from the USA. The number-average molecular weight (Mn) is an important reference value for the construction of the average molecular structure of residue oil, representing the contribution of each component to the system. The formula for Mn is as follows: ; Where H i It is relative to the retention volume V i The baseline detector signal level; M i It is the retained volume V i The molecular weight of the polymer component at the given location; n is the number of data points.
[0025] The number-average molecular weight (Mn) of vacuum residue was measured to be 826 g / mol. Then, 826 g / mol will be used as the reference molecular weight to calculate the approximate number of each atom, and corrections will be made based on the performance of the residue.
[0026] Elemental Analysis: The elemental analyzer used was the Elementar Isoprime PrecisION. Test details are as follows: 1. Test environment: temperature 22℃, humidity 20%. 2. Test modes: CHNS mode and O mode. 3. Equipment test conditions: ashing tube temperature 1150℃, reduction tube temperature 830℃, pyrolysis tube temperature 1170℃. 4. Test quality control method: The accuracy of the test results was monitored using four tests with standard substances, and the precision was monitored by 100% repeatability of the sample. 5. Quality control evaluation: The relative deviation between two repeat measurements was less than 2%. The elemental content of the vacuum residue is shown in Table 1 as a mass fraction.
[0027] .
[0028] Atom count calculation: Determine the mass fraction of each element, and calculate the approximate number of each element using the mass fraction and the above Mn=826g / mol: The above rounded results are shown in Table 2: ; Functional groups were determined using Fourier Transform Infrared Spectroscopy (FT-IR) with an instrument model INVENIO R. Liquid testing of the residue oil was performed, and the spectral data of the residue oil sample obtained simultaneously are as follows: Figure 1 As shown.
[0029] according to Figure 1 It can be seen that the characteristic functional groups of vacuum residue appear at wavenumbers of 3000~650cm. -1Nine distinct characteristic peaks were observed; the wavenumber positions of these characteristic peaks were 2920, 2850, 1603, 1460, 1375, 867, 812, and 712 cm⁻¹. -1 Table 3 provides an analysis of the characteristic peaks for each wavenumber: .
[0030] Therefore, according to Figure 1 Table 3 shows its chemical structure as follows: 1) Carbon exists mainly in the form of aromatic structural framework and methylene form, with aromatic ring skeleton as the main form, followed by side chains that replace aromatic structures, and a small amount in the form of methyl.
[0031] 2) The oxygen functional group mainly exists in the form of furan aromatic heterocycles.
[0032] 3) Sulfur exists mainly in the form of thiophene aromatic heterocycles in the functional groups.
[0033] 4) Quantitative analysis of molecular structure shows that vacuum residue has a high content of alkyl side chains (-CH2-, -CH3), residual aliphatic structure, and cycloalkane structure.
[0034] Furthermore, based on the Mn concentration of 826 g / mol obtained by gel permeation chromatography, and after appropriate rounding of the data obtained by Fourier transform infrared spectroscopy structural analysis, the molecular formula of vacuum residue was finally obtained as C. 57 H 74 OS, as shown in Table 4.
[0035] .
[0036] Step S2: Construct the average molecular structure model of vacuum residue, including the following steps: Step S2.1, Core Structural Principles: Based on experimental infrared spectroscopy analysis and the low H / C ratio, it can be seen that the structure contains benzene rings and cycloalkanes. The vast majority of cycloalkanes in the residue oil are cyclohexanes, thus forming a structure dominated by benzene rings and cyclohexane. Since the H of the benzene ring is hardly reflected in the experimental infrared spectrum, it can be seen that the H of the benzene ring is largely substituted, and the benzene ring preferentially links to the side chain.
[0037] Step S2.2, Functional Groups and Side Chain Hypotheses: 1) Assume that each cycloalkane ring and benzene ring is connected to an alkyl side chain, and there are no multiple branches.
[0038] 2) The oxygen functional group exists in the form of a furan aromatic heterocycle. Since there is only one oxygen group, there is only one furan group.
[0039] 3) Sulfur exists mainly in the form of thiophene aromatic heterocycles in the functional groups.
[0040] Step S2.3, Construction of the vacuum residue model: First, the structural derivation is performed, which includes the following: Calculate the degree of unsaturation: ; Where: C: number of carbon atoms; H: number of hydrogen atoms; N: number of nitrogen atoms (0 here); X: Number of halogen atoms (0 here); Oxygen (O), Sulfur (S): Do not affect the calculation of unsaturation, no additional processing is required in the formula; Ω=21, which limits the upper limit of the total number of rings and double bonds in the molecule, providing constraints for subsequent structure allocation.
[0041] Next, structural features are derived, including the following: Based on the qualitative evidence of the aromatic structure by infrared spectroscopy, and following the typical structural characteristics of heavy oil with "polycyclic fused aromatics as the core", assuming that the fused aromatic ring and the heterocycle do not share double bonds, the ring structure was verified by multiple schemes of unsaturation matching, and the results are shown in Table 5.
[0042] .
[0043] Among the above ring number combinations, the scheme with 4 fused aromatic rings, 3 cyclohexyl rings, and 2 heterocycles is highly consistent with the total unsaturation degree Ω=21 and conforms to the typical structural characteristics of "partially hydrogenated aromatics" in heavy oils. Considering that the fused aromatic rings and heterocycles inevitably share conjugated double bonds in this system, which would lead to double counting of the unsaturation degree, and that this fusion form is more common in heavy oil molecules, the ring number distribution was modified, and the final structural scheme was determined to be 4 fused aromatic rings, 5 cyclohexyl rings, and 2 heterocycles.
[0044] Next, the length and position of the side chains are determined, specifically including the following: Long-chain alkyl groups (2920, 2850 cm⁻¹) in infrared spectra -1 The strong absorption indicates the presence of long straight-chain side chains. Given a total carbon number of 57, after deducting the core fused ring (approximately 2 carbons), the remaining 28 carbons are distributed into 3 long side chains. These side chains are attached to the outside of the aromatic / cycloalkane ring, thus avoiding disruption of the aromatic conjugated system.
[0045] Finally, molecular model construction and visualization are performed, which includes the following: Three-dimensional molecular modeling was completed using Materials Studio, and two-dimensional structural formulas were drawn using ChemSketch. The resulting structural model is shown below. Figure 2 As shown.
[0046] Step S3: Generate 3D structure file: including the following steps: First, by exporting the file and converting it to an ORCA .inp file, write the most basic ORCA-compatible .inp file as follows: First, fill in 1. DFT method; 2. Dispersion; 3. Diffuse; 4. Basis set; 5. Run content; then fill in the parallel method (modify the script information if calculating on a platform); next, fill in the implicit solvent (not used in this article); since this system is neutral, the total charge in the file must be 0. At this time, each atom either does not fill in charge information or the sum of charges is 0, otherwise an error will occur; then fill in the DFT spin multiplicity. When all molecular orbitals are filled with paired electrons, the electronic structure system without empty single-occupied orbitals or unpaired electrons is also called a closed-shell system. At this time, the number of unpaired electrons n=0, and the multiplicity is 1 (singlet state); when a single free radical exists, the number of unpaired electrons n=1, and the multiplicity is 2. Finally, provide the atomic coordinate information to complete the ORCA .inp file, and finally optimize the configuration ball-and-stick model as follows. Figure 3 As shown: Then, the Optimized Potential Energy (OPLS-AA) force field of the all-atom version liquid simulation was added using Autoff software, because the force field basically matches the experimental data in the main properties and temperature relationship of the asphalt mixture model, which provides a basis for vacuum residue oil that is very similar to asphalt. QEQ charge was assigned, and moltemplate software was selected for output to obtain pdb and lt files.
[0047] Secondly, create a packmol.inp file, and then use packmol to randomly generate 20 molecules within a certain range by inputting the pdb file mentioned above, to obtain a new pdb file. At this time, although the molecules are randomly arranged, they do not overlap.
[0048] Finally, use moltemplate to assign boundaries and output a complete data file. Take the lt file obtained above as the input file, set the box size to 50Å×50Å×50Å, and finally obtain a data file with OPLS-AA force field information. Copy the other files in the order required in the run.in.EXAMPLE file (Init Section, Atom Definition Section, Settings Section). Then add the ensemble to be run to obtain the complete in file.
[0049] II. Validation of the vacuum residue model: Step S4: Density Verification: NPT ensemble simulations are performed in Lammps. Lammps is a mainstream dynamics simulation software that includes numerous force fields, such as the classical all-atomic force field CHARMM, the OPLS-AA force field, the metallic force field EAM, and the reactive force field ReaxFF. The "in" and "data" fields mentioned above are for the subsequent use of Lammps.
[0050] First, the energy of the constructed residue oil system is minimized using the conjugate gradient method (Cg), with the convergence criterion being that the interatomic force is less than 10. -8 The system energy change is less than 10. -10 The maximum number of iteration steps and the maximum number of force calculation steps are both set to 100,000. The purpose is to eliminate unreasonable high-energy stress, atomic overlap, bond length and bond angle distortion in the initial structure, and reduce the system to a stable configuration with the lowest local energy.
[0051] Secondly, the system is annealed under an NTP ensemble (isothermal-isobaric ensemble): the system is heated from 300K to 600K at a constant heating rate, held at 600K for a sufficient period of isothermal activity to activate molecular chain segment motion, and then slowly cooled to 300K. Annealing can release the large stresses caused by long-range molecular entanglement, interlayer compression, and local distortion of O / S heteroatoms.
[0052] Then, to eliminate porosity and molecular packing defects generated during the modeling process, the system was subjected to pressure cycling at 300 K under the NPT ensemble: the pressure was gradually increased from 1 bar to the target pressure of 1000 bar, and equilibration was maintained under high pressure to compress system voids and optimize molecular packing, followed by gradual decompression back to 1 bar. This operation significantly improved the density uniformity of the system, making the density of the simulated system closer to the actual experimental value, and obtaining a more thermodynamically stable dense structure.
[0053] Finally, a longer NP ensemble simulation was performed at 300K and 1 bar, such as... Figure 4 The image shows the density distribution of vacuum residue in the last 50 ps, with an average value of 0.99554 g / cm³. 3 The experiment was conducted according to GB1884 standard, and the average of two results was taken: 0.99405 g / cm³. 3 Error: 0.1499%, which is within the acceptable range for molecular dynamics simulations.
[0054] Step S5: Infrared spectroscopy verification.
[0055] Theoretically calculated infrared vibrational wavenumbers are typically higher than experimentally measured values due to differences in functional approximation, basis set truncation, and quantum mechanical descriptions of vibrational modes. To quantitatively match the calculated spectra with experimental results, this embodiment uses the same functional, basis set, and dispersion conditions (B3LYP D3BJ def2-SVP) to measure C5H...12 C 10 H 22 C 15 H 32 C 20 H 42 Wavenumber correction was performed, but the relative intensity of the absorption peaks was not adjusted. For example... Figure 5 The image shows the effect of C5H 12 C 10 H 22 C 15 H 32 C 20 H 42 The infrared spectrum used for calculation, such as Figure 6 The image shown is a magnified view of the high wavenumber region, as follows: Figure 7 The image shown is a magnified view of the low wavenumber region. The conclusions are as follows: 1. Methyl antisymmetric stretching vibration at wavenumber 3097 cm⁻¹ -1 Nearby, corresponding to the experimental wavenumber of 2950 cm⁻¹ -1 The correction factor is 0.9525; the methylene antisymmetric stretching vibration at wavenumber 3065 cm⁻¹ -1 Nearby, corresponding to the experimental wavenumber of 2920 cm⁻¹ -1 The correction factor is 0.9526; the correction factor for antisymmetric stretching vibration in the high-frequency region is 0.9525.
[0056] 2. The symmetric stretching vibrations of methyl and methylene groups at wavenumber 3016 cm⁻¹ -1 Nearby, the peaks in the DFT calculation overlapped, and the intensity increased, corresponding to experimental wavenumbers of 2870 and 2850 cm⁻¹. -1 The average value is 2860cm. -1 The correction factor is 0.9483, therefore the correction factor for symmetrical stretching vibration in the high-frequency region is 0.9483.
[0057] 3. Methylene scissor bending vibration and methyl asymmetric bending vibration at wavenumber 1471 cm⁻¹ -1 Nearby, corresponding to the experimental wavenumber of 1450 cm⁻¹ -1 The correction factor is 0.9857; the methyl symmetric bending vibration at a wavenumber of 1404 cm⁻¹ -1 Nearby, corresponding to the experimental wavenumber of 1375 cm⁻¹ -1 The correction factor is 0.9793; the cooperative vibration of (CH2)n (n>4) at wavenumber 740 cm⁻¹ -1 Nearby, corresponding to the experimental wavenumber of 720cm. -1 The correction factor is 0.9729; the overall correction factor for vibration in the low-frequency region is 0.9793.
[0058] The DFT frequency calculation first checks for the presence of imaginary frequencies; if imaginary frequencies are found, the suitability of the pseudopotential needs to be checked. The infrared (IR) absorption spectrum is obtained through the DFT frequency calculation and the aforementioned correction factors. This helps verify the reliability of the constructed vacuum residue model by comparing the calculated IR absorption spectrum with the aforementioned... Figure 1 The spectra obtained from vacuum residue experiments were compared, such as... Figure 8 As shown.
[0059] As shown in the experimental infrared spectra, the calculated infrared spectra are consistent with the experimental data. Regarding peak positions, the characteristic absorption bands predicted by DFT match well with the experimental wavenumbers, reflecting the accurate reproduction of molecular vibrational modes in the vacuum residue model. However, regarding peak intensity, due to the lack of intensity correction, the relative intensities of the calculated peaks also show some differences from the experimental observations.
[0060] III. Pyrolysis Reaction Simulation: Step S6: Construction and verification of the pyrolysis reaction model.
[0061] First, a structural model of vacuum residue was established by combining experimental characterization and theoretical calculation. The pyrolysis process was simulated using ReaxFF-MD with the force field parameter CHONSSi.ff, which is suitable for simulating vacuum residue pyrolysis. The model was relaxed for 50 ps under the NVT ensemble at 300 K using a Nose-Hoover thermowell with a damping time constant of 100 fs and a step size of 0.10 fs. The lattice edge length parameters are shown in Table 6. Using a smaller density saves computational effort, avoids atomic loss, and has minimal impact on the reaction mechanism.
[0062] .
[0063] Force field parameter verification: ReaxFF compares the force field, charge, and DFT's mayer bond order and CM5 charge, as shown below. Figure 9 , 10 As shown: the typical bond order of a single atom ranges from 0 to 3, with overall values on the order of 0.1. Using the lowest bond order H atom (0.9) as a reference, Figure 9 The average absolute error (MAE) of the bond order comparison is 0.0463, which is at a low level. The local single-point deviation is small, and the quantitative accuracy of the force field is excellent. R2=0.814 indicates that the linearity of the bond order simulated by ReaxFF and the bond order calculated by DFT is at a medium to high level, which can support the system simulation calculation. The fitting slope is 0.857, and the overall fitting line is slightly lower than that of the DFT-mayer bond order.
[0064] Figure 10The mean absolute error (MAE) of the charge comparison is 0.0166e, which is at an extremely low level, with very small quantitative deviation. The numerical prediction accuracy of the force field for atomic charge is excellent. The R² = 0.892 is close to 0.9, indicating excellent linear correlation. This shows that the charge distribution calculated by ReaxFF-QEQ is highly consistent with the overall charge change trend calculated by DFT-CM5 quantization, and the force field can well reproduce the charge distribution law of the system. The fitting slope is 1.102, and the overall fitting line is slightly higher than that of DFT-CM5 charge.
[0065] The cross-validation of the two core microscopic parameters, bond order and charge, shows that the bond order description capability of the CHONSSi.ff force field and the charge calculation accuracy of ReaxFF-QEQ both exhibit good consistency with the DFT results, fully demonstrating the reliability of the force field parameters and the applicability of the system, and enabling it to be applied to the reaction molecular dynamics simulation study of the vacuum residue model.
[0066] Step S7: Analysis of pyrolysis mechanism.
[0067] This application's embodiments explore the structural evolution characteristics of a vacuum residue model under different heating rates. The simulation uses a time step of 0.1 fs. While this step increases computation time, it accurately captures reaction path details and effectively reduces the probability of atom loss under high-temperature conditions. In Reaction Force Field (ReaxFF) simulations, atom loss is a common cause of computational interruption; by default, atom loss directly triggers an error and terminates the simulation. To ensure simulation continuity, it is recommended to configure the `in` file to periodically save a restart file, for example, setting a restart file to be generated every 20,000 steps (restart20000tmp.restart). If atom loss occurs, the most recent restart file can be loaded, and the time step can be reduced to rerun the simulation, thus reducing the risk of further atom loss. NVT equilibration was performed using a force field of 300K, 30ps, and CHONSSi.ff. The system spontaneously evolved to a stable, low-stress equilibrium state that met the set temperature and pressure requirements, while maintaining the integrity of the vacuum residue molecular structure. Energy minimization was achieved using the conjugate gradient method (cg). cg utilizes historical gradient information to construct conjugate directions, avoiding repeated searches and resulting in faster convergence. The energy variation tolerance and maximum force component tolerance values are both 10. -8 For the 20 large molecule model, intermolecular packing optimization and energy convergence can be achieved in 30 ps. Both relaxation and production simulations used an NVT ensemble with a Nose-Hoover temperature controller (damping constant: 100 fs). This method strictly satisfies microcanonical ensemble statistics and is suitable for long-term equilibrium and production simulations.
[0068] The details are explained below.
[0069] Step S7.1, Reactant Data Analysis: To reveal the regulatory effect of heating rate on the pyrolysis kinetics of vacuum residue oil, and to obtain the product yields at different temperatures, the final number of steps was determined by the difference between the final and initial temperatures and the heating rate. Heating was performed at three rates: 10 K / ps, 12 K / ps, and 16 K / ps, with the final temperature kept constant at 3500 K. The following figure shows the yields of each component at different heating rates: Figures 11, 12, and 13 show the curves illustrating the changes in the proportions of each component in vacuum residue at different heating rates. The pyrolysis behavior at different heating rates indicates that the reaction process exhibits significant temperature-rate coupling characteristics and can be divided into four stages: 1) First stage (300–2000 K): The reactants mainly undergo side chain breakage. At this time, the molecular thermal motion energy level is low, and the breaking of macromolecular bonds requires overcoming a high activation energy barrier, resulting in a lag in the pyrolysis reaction.
[0070] 2) Second stage (2000–2500K): As the heating rate increases, the frequency of molecular collisions increases, which effectively reduces the activation energy of the reaction, accelerates the pyrolysis rate, generates a large number of gaseous products, and the residue oil begins to decompose into small molecules.
[0071] 3) The third stage (2500–3000K): The pyrolysis rate reaches its peak, and the gaseous products are mainly generated in this stage; the residual oil gradually completes the reaction and condenses into coke macromolecules, the proportion of heavy oil increases, part of the light oil is converted into gas, and the other part condenses with the heavy oil into coke.
[0072] 4) Fourth stage (3000–3500K): The proportion of heavy oil and light oil decreases, gas content increases slightly, and carbon content rises. 40+ The coke tends to stabilize, and the system enters the pyrolysis rate saturation zone.
[0073] In terms of component proportions, the reaction trends at different heating rates are generally consistent: at 10 K / ps, due to the longer reaction time, the formation of light products is earlier, and the condensation reaction is more complete after 3000 K, the proportion of coke does not decrease, and the reaction mechanism is slightly different; at 12 K / ps and 16 K / ps, the initial reaction time is shorter, and there are still more residual oils at 2000 K and 2500 K, and the product composition at different heating rates tends to be similar after 3000 K.
[0074] Thermogravimetric behavior analysis: Figure 14 The relationship between temperature and residual reactant mass in the actual pyrolysis and molecular simulation of vacuum residue was presented. Under experimental conditions, vacuum residue showed almost no pyrolysis reaction from room temperature to 340℃, with the main thermal weight loss concentrated in the range of 340℃ to 502℃. During this stage, macromolecular decomposition was intense, and the reactant mass decreased rapidly. When the temperature exceeded 502℃, the pyrolysis rate slowed down significantly, entering a later stage where deep decomposition and condensation polymerization coexisted.
[0075] In molecular simulation studies, to overcome the bottleneck of extremely slow reaction kinetics and low computational efficiency at low temperatures, a high-temperature acceleration strategy was adopted: the simulation started by triggering the reaction at 1200K, with the main pyrolysis range from 1200K to 2500K. Although the simulation temperature is much higher than the experimental temperature, increasing the temperature in the kinetic simulation does not change the core reaction mechanism, thus significantly shortening the simulation time and saving computational resources while ensuring consistency in physical meaning.
[0076] Overall, the experimental and simulated pyrolysis curves show a high degree of agreement, both exhibiting a three-stage characteristic of "inert heating – rapid pyrolysis – slow reaction." The slight unevenness in the simulated curves can be attributed to the cutoff radius setting in the molecular simulation: even if chemical bonds are broken, if the product molecules do not exceed the cutoff radius, they will still be considered unreacted by the system, leading to slight fluctuations in the mass change curve. This does not affect the judgment of the overall pyrolysis trend and reaction mechanism.
[0077] In summary, molecular simulations conducted using a high-temperature acceleration strategy can effectively reproduce the macroscopic thermogravimetric behavior of vacuum residue pyrolysis, providing a reliable computational basis for in-depth analysis of its microscopic reaction mechanism.
[0078] Key product analysis: The generation data of ethylene and hydrogen, the beneficial products of vacuum residue, were compared and analyzed at different heating rates and final temperatures.
[0079] Ethylene (C2H4): Figure 15 The effects of different heating rates and pyrolysis temperatures on the C2H4 product distribution in a vacuum residue model were demonstrated. Data showed that after 2070 K, the number of ethylene molecules gradually increased with increasing pyrolysis temperature, reaching peak values at 2964 K (10 K / ps), 3234 K (12 K / ps), and 3274 K (16 K / ps), respectively. Subsequently, with further heating, the number of ethylene molecules decreased, decreasing by approximately 20 molecules by 3500 K. Based on this, it can be inferred that the optimal simulation temperature range for vacuum residue shifts with increasing heating rate, remaining at a relatively high level between 2900 K and 3300 K.
[0080] hydrogen: Figure 16 The effects of different heating rates and pyrolysis temperatures on the H2 product distribution in a vacuum residue model were demonstrated. The results showed that after 2289 K, the hydrogen production gradually increased with increasing pyrolysis temperature, and did not reach its peak before 3500 K; the hydrogen production trend was highly consistent under different heating rates. This is because hydrogen is generated by the combination of two hydrogen free radicals (2H·→H2), and a large number of hydrogen free radicals are continuously generated during pyrolysis, with a uniform generation mechanism, thus resulting in a continuous upward trend in hydrogen production.
[0081] Heteroatom transformation pathways: such as Figure 17 and Figure 18 As shown, at lower temperatures (1500K–2500K), the main reaction in vacuum residue is the scission of alkyl side chains. The broken side chains continue to break down, forming smaller free radicals and stable substances. At higher temperatures (above 2500K), unreacted vacuum residue molecules have a chance to directly open their rings or split off smaller transition substances from their side chains. The large amount of H· produced during the reaction will react with other free radicals to form more stable substances, including reacting with H· itself to form H2.
[0082] The pathway of ethylene production from vacuum residue via pyrolysis reveals that at lower pyrolysis temperatures, the alkyl side chains of polycyclic aromatic hydrocarbons (PAHs) and cycloalkanes react first, decomposing to form C2H5·. C2H5· then reacts with H· to generate C2H4. At higher pyrolysis temperatures, some carbon-carbon bonds on the rings are broken, while new unsaturated bonds are formed. The intermediates further undergo pyrolysis, with the side chains and carbon-carbon bonds on the rings continuing to break, resulting in the release of smaller molecular fragments and ultimately the formation of ethylene (C2H4). This process involves the gradual breaking, dehydrogenation, and miniaturization of PAHs and cycloalkanes, promoting the breaking and recombination of bonds on the rings.
[0083] From the reaction mechanism of O (and S similarly), it can be seen that for O to become H2O, HO· and other substances, the furan aromatic heterocycle needs to open its ring to become an alcohol. O eventually produces H2O and alcohols, while S produces very little H2S. This is one of the reasons why high-sulfur residue oil has poor quality.
[0084] Temperature has a significant regulatory effect on the pyrolysis reaction pathway of vacuum residue. In the low-temperature range of 1500–2500 K, vacuum residue mainly undergoes alkyl side chain scission reaction, and the alkyl free radicals generated by scission will continue to decompose, generating small molecule free radicals and stable products. When the temperature rises to the high-temperature range above 2500 K, unreacted large molecules of vacuum residue can directly undergo ring-opening of the cyclic skeleton, or generate small molecule transition intermediates through side chain decomposition; the large number of hydrogen free radicals (H·) generated during the reaction can combine with various free radicals to form stable compounds, and at the same time, hydrogen free radicals couple with each other to generate hydrogen gas (H2).
[0085] Analysis of the ethylene formation pathway reveals that during the low-temperature pyrolysis stage, the alkyl side chains of polycyclic aromatic hydrocarbons (PAHs) and cycloalkanes preferentially break, generating ethyl radicals (C2H5·). These ethyl radicals further combine with hydrogen radicals to form ethylene (C2H4). During the high-temperature pyrolysis stage, the carbon-carbon bonds in the cyclic skeleton break and form unsaturated bonds. The intermediates continue to cleave, and small molecular fragments are continuously broken and shed from the side chains and intracyclic carbon-carbon bonds, ultimately forming ethylene. Overall, the ethylene formation process involves the skeletal breaking, dehydrogenation, and miniaturization of PAHs and cycloalkanes, achieving the breaking and recombination of cyclic chemical bonds.
[0086] Regarding heteroatom reactions, oxygen and sulfur exhibit similar reaction patterns. Oxygen-containing compounds require ring-opening of furan aromatic heterocycles to generate alcohol intermediates, which are then converted into products such as H2O and hydroxyl radicals (HO·). Ultimately, oxygen mainly exists in the form of water and alcohols. Sulfur-containing thiophene heterocycles are more stable, generating only trace amounts of hydrogen sulfide (H2S) during pyrolysis. These trace amounts of harmful sulfur-containing gases are also one of the important reasons for the poor quality of high-sulfur vacuum residue oil.
[0087] In this embodiment, it should also be noted that the method further includes: IV. Analyzing the structural characteristics of vacuum residue based on the results of molecular dynamics simulation of pyrolysis reaction.
[0088] Vacuum residue differs from light oils in molecular composition and microstructure, and its structural characteristics directly determine its pyrolysis reactivity. This example study determined the molecular structure of the residue to be C2. 57 H 74 The OS sample has a carbon content of 82.9% by mass, characterized by high carbon content, abundant alkyl side chains (-CH2-, -CH3), and a high proportion of aliphatic structure residues. The heteroatoms inside the sample are present in the form of stable aromatic heterocycles, with oxygen present in furan heterocycles and sulfur present in thiophene heterocycles; the molecular skeleton is mainly composed of polycyclic aromatic hydrocarbons and cyclohexane structures, and generally contains side chain substituents.
[0089] This embodiment primarily combines DFT, ReaxFF molecular dynamics (ReaxFF-MD) simulations, and thermogravimetric experiments to explore the pyrolysis behavior of vacuum residue under different operating conditions. By optimizing modeling methods and simulation schemes, the macromolecular modeling strategy and key simulation parameters are clarified, revealing the formation rules of pyrolysis products and elucidating the microscopic mechanisms of product formation. At the modeling level, an accurate macromolecular model of vacuum residue is constructed, and optimal simulation parameters are determined, filling the research gap in accurate molecular models of high-sulfur vacuum residue and improving the model's adaptability and practicality in industrial simulation. At the reaction mechanism level, the structure-activity relationship study of high-sulfur vacuum residue is improved, deepening the understanding of the mechanisms of heteroatom transformation and macromolecular fragmentation. Furthermore, the above can provide a theoretical basis for the targeted regulation of the vacuum residue pyrolysis process and provide technical support for promoting the resource utilization of vacuum residue and expanding alternative oil resources.
[0090] It should be noted that the above-described embodiments are merely illustrative of several implementation methods of this application, and while the descriptions are relatively specific and detailed, they should not be construed as limiting the scope of the patent application. It should be pointed out that, for those skilled in the art, other equivalent embodiments or equivalent technical means may be included based on the concept and disclosed technical solutions of this application, and all such embodiments or means fall within the protection scope of the claims of this application.
Claims
1. A method for constructing and validating a vacuum residue model, characterized in that, Includes the following steps: Step 1: Obtain experimental data of vacuum residue, including elemental analysis data, number-average molecular weight data, and infrared spectral data; Step 2: Construct an average molecular structure model of vacuum residue based on the experimental data; Step 3: Perform density verification of the average molecular structure model using molecular dynamics simulations, including performing energy minimization, annealing, and pressure cycling on the system in sequence, and calculating its density under isothermal and isobaric ensembles, and comparing the calculated results with the experimental density for verification. Step 4: Perform quantum chemical calculations to verify the constructed and validated vacuum residue model, including calculating its infrared spectrum, Mayer bond order and CM5 charge, and compare and verify the results with experimental infrared spectra and reaction force field simulations. Step 5: Based on the validated vacuum residue model, construct a pyrolysis reaction model and use reaction force field molecular dynamics simulation to study its pyrolysis behavior.
2. The method according to claim 1, characterized in that, Step 2, which involves constructing the average molecular structure model of vacuum residue, specifically includes: Based on elemental analysis data and number-average molecular weight data, the molecular formula of vacuum residue was determined to be C2. 57 H 74 OS; Based on the calculation of unsaturation degree and infrared spectral analysis, the molecular structure was determined to be based on a core skeleton of fused aromatic rings, cyclohexane rings, oxygen-containing heterocycles, and sulfur-containing heterocycles, with alkyl side chains attached. Based on the assumption that each cycloalkane ring and benzene ring is connected to an alkyl side chain and there are no multiple branches, two-dimensional and three-dimensional average molecular structure models are constructed.
3. The method according to claim 1, characterized in that, The density verification in step 3 specifically involves: The energy of the system was minimized using the conjugate gradient method. Then, heating-cooling annealing and pressure cycling were performed sequentially under isothermal and isobaric ensemble conditions. Finally, isothermal and isobaric ensemble simulations were performed at 300 K and 1 bar to calculate the system density and compare it with the experimental density. If the relative error was within the acceptable range of molecular dynamics simulation, the model passed density verification.
4. The method according to claim 1, characterized in that, The quantum chemical calculation verification in step 4 specifically includes: The vibrational frequencies of the constructed model were calculated using density functional theory, and the calculated wavenumber was corrected. The peak positions and peak intensities of the corrected calculated infrared spectrum were compared with those of the experimental infrared spectrum. The Mayer bond order and CM5 charge of the calculated model were compared with the bond order and charge obtained from the reaction force field simulation. The mean absolute error and coefficient of determination were calculated to verify the reliability of the force field parameters.
5. The method according to claim 4, characterized in that, The specific steps for correcting the calculated wavenumber include: A correction factor of 0.9525 is used for antisymmetric stretching vibration in the high-frequency region, a correction factor of 0.9483 is used for symmetric stretching vibration in the high-frequency region, and a correction factor of 0.9793 is used for vibration in the low-frequency region.
6. The method according to claim 1, characterized in that, Step 5 specifically includes: Based on the validated vacuum residue model, a periodic pyrolysis reaction model containing multiple molecules was constructed. Reaction force field molecular dynamics simulations were performed at different heating rates to analyze the residual mass of reactants, the formation rules of products, and the evolution trend of key products ethylene and hydrogen during pyrolysis. Based on the simulation results, the pyrolysis process is divided into three stages: inert heating, rapid pyrolysis, and slow reaction. The effects of different heating rates on product distribution are analyzed.
7. The method according to claim 6, characterized in that, The different heating rates include 10K / ps, 12K / ps and 16K / ps, with a final heating temperature of 3500K. The optimal simulated temperature range for vacuum residue pyrolysis is determined to be delayed as the heating rate increases, with ethylene yield remaining at a relatively high level in the range of 2900K to 3300K.
8. The method according to claim 1, characterized in that, It also includes step 6: Based on the results of molecular dynamics simulation of pyrolysis reaction, the structural characteristics of vacuum residue were analyzed. These structural characteristics include: high carbon content, abundant alkyl side chains, high proportion of aliphatic structure residues, heteroatoms in the form of aromatic heterocycles in furan and thiophene rings, and a molecular skeleton dominated by polycyclic aromatic hydrocarbons and cyclohexane structures.
9. A system for constructing and validating a vacuum residue model, characterized in that, include: The data acquisition module is used to acquire experimental data of vacuum residue, including elemental analysis data, number-average molecular weight data, and infrared spectral data. The model building module is used to construct an average molecular structure model of vacuum residue based on the experimental data. The density verification module is used to verify the density of the average molecular structure model through molecular dynamics simulation. The quantum chemical verification module is used to verify the average molecular structure model by infrared spectroscopy, Mayer bond order and CM5 charge through quantum chemical calculations. The pyrolysis simulation module is used to build pyrolysis reaction models based on validated models and perform molecular dynamics simulations of reaction force fields.