Unsteady-State Mathematical Modeling of Hydrocarbon Feedstock Pyrolysis

: Hydrocarbon feedstock pyrolysis is an important method for obtaining monomers that are then used to produce various polymer materials. During this process, a mixture of hydrocarbons is heated at a high temperature and in the absence of oxygen. Because of the side reactions of polymerization and polycondensation, coke products are formed and settle on the inner walls of the coil. This decreases the technical efficiency of the hydrocarbon pyrolysis furnace during its operation, making the process unsteady. In the present research, we developed an unsteady-state mathematical model of hydrocarbon feedstock pyrolysis in order to improve the monitoring, forecasting, and optimization of this technological process. This model can calculate the rate of coke deposition along the length of the coil, considering the technological parameters and the composition of the supplied raw materials (the calculated value of coke deposition rate equals 0.01 mm/day). It was shown that with an increase in the propane/butane ratio from 4/1 to 1/4 mol/mol, the ethylene concentration decreases from 3.45 mol/L to 3.35 mol/L.


Introduction
Currently, almost all the objects that surround people and are used in everyday and professional activities are made from various polymer materials. Among these, the most common materials include polyethylene, polypropylene, polyvinyl chloride, polyethylene terephthalate, polystyrene, and many others. The substances used to produce these and many other polymers are the alkene monomers ethylene and propylene. Each of these alkene monomers has a single double bond that is responsible for their activity in polymerization reactions [1,2].
Thus, the high and unabated demand in the petroleum product market for ethylene and propylene obviously imposes rather stringent requirements not only on the quality of semi-finished products, but also, more importantly, on their quantity. This largely depends on the composition of the processed raw material, which, in turn, determines the optimal modes of the pyrolysis stage, which is the main method for obtaining the necessary monomers. During this process, a mixture of hydrocarbons is heated at a high temperature and in the absence of oxygen.
Today, one of the most common designs for the pyrolysis process is a tubular furnace. Devices of this type usually include two sections. The first one, located at the top of the apparatus, is convection. Here, the feedstock supplied to the apparatus is heated by the heat of the flue gases, evaporates, and is mixed with water vapor, also supplied to the apparatus, while being heating to the temperature at which the decomposition reaction will occur. The supplying of water vapor is necessary, since this makes it possible to reduce the rates of side reactions during pyrolysis. As a result, the partial pressure of hydrocarbons decreases and this facilitates the occurrence of reactions that increase the volume, that is, targeted primary reactions of hydrocarbons decomposition. The second, lower section is radiant. In this section, the heat generated by combustion in the fuel burners is already needed only to support the occurrence of hydrocarbon decomposition reactions. It should be noted that there are from two to eight coils inside the furnace in which these reactions take place.
The technical efficiency of the hydrocarbon pyrolysis furnace is determined both by the composition of the feedstock and by the thermodynamic parameters of the process. However, it inevitably decreases during the operation of the installation. This is due to formation of coke deposits during the pyrolysis process. These deposits are the result of the pyrolysis side reactions of polymerization and polycondensation, the products of which settle on the inner walls of the coil. This decreases the section's cross section and, in turn, increases in the pressure drop at the ends of the tube. As a result, the contact time increases, and the coke layer thickness becomes even thicker. In addition, coke deposits on the coil walls can create areas of thermal stress. This can lead to burnout, which is already an emergency situation that may require replacing of the entire coil.
According to the stated problems, hydrocarbons pyrolysis modeling is worth optimizing. There are two main directions of improving the high-temperature process of hydrocarbon feedstock pyrolysis. First is the modernization of the technology itself [3][4][5][6][7][8][9]. Within the framework of this direction, the burner devices are being modernized [10][11][12][13], research is being carried out to initiate the pyrolysis process, new and more efficient catalysts are being developed. The second direction is using the mathematical modeling method to search for optimal process modes. Both deterministic and stochastic models are developed and used [14][15][16][17]. A common feature of the problems solved in modeling the pyrolysis process is the constancy of the feed stock composition. However, in real conditions, especially in the case of supplying pyrolysis feedstock simultaneously from several enterprises, there are strong variations in the feedstock composition.
At present, there are a large number of studies of the pyrolysis process itself: the study of chemistry, kinetics, thermodynamic and hydrodynamic laws, and the development with further improvement of its mathematical models. Nevertheless, in these works, the processes of coke formation, its precipitation, and combustion are not considered. These processes narrow the performed developments and complicate the forecasting of the technological process. Having significant experience in unsteady processes modeling [18][19][20][21], we developed an unsteady-state mathematical model of hydrocarbon feedstock pyrolysis in order to improve monitoring, forecasting, and optimization of this technological process. As a base, we used the analysis of the hydrocarbon pyrolysis process and a previously proposed kinetics model, which is based on a somewhat simplified scheme of the process, allowing, nevertheless, for the calculation of concentrations of the main components of the reaction mixture for a wide range of process temperatures [17,22].
There are at least three main and significant problems that arise during development of mathematical models of technological processes. First, the volume of necessary calculations and processed information grows quickly, since determination of a large number of numerical values of the physicochemical parameters of a model is a nontrivial task. Second, the analysis of the results of calculations obtained using a mathematical model can be to a certain extent complicated by the large number of components of the reaction mixture, since, as a rule, the models do not consider their additional effects of interaction with each other. Thirdly, it is very difficult to compare the calculated data obtained using the mathematical model and the experimental data obtained from the existing production.
According to the information above, the aim of the study was to predict the operation of the pyrolysis unit, considering the changes in technological parameters, the composition of raw materials, and the amount of deposited coke.
To achieve this goal, it was necessary to perform a number of tasks: (1) study the mechanism of the pyrolysis process itself and the mechanism of coke deposition; (2) perform a thermodynamic analysis of the pyrolysis process; (3) develop a formalized reaction network according to the known mechanism considering the thermodynamic analysis of the pyrolysis process; (4) develop a kinetic model of the pyrolysis process; (5) develop an unsteady mathematical model of a pyrolysis reaction furnace with due consideration of the coke deposition on the inner surface of the coil; (6) complement the system of equations with differential equations for changing the thickness of the coke layer and the pressure of the reaction medium along the length of the coil with the astronomical time; (7) solve the inverse kinetic problem; (8) verify the mathematical model; (9) perform predictive calculations for the duration of the cycle of the pyrolysis furnace; (10) perform optimization calculations to obtaining the maximum degree of conversion, selectivity, and duration of the cycle of the pyrolysis reaction furnace.
The novelty of this research lies in developing an unsteady mathematical model of the pyrolysis process. This model considers deposition of coke on the inner surface of the coil. The rate of coke deposition depends on the feedstock composition, on the temperature and pressure parameters, and on the mode of the medium movement in the apparatus. Such a model for the pyrolysis process has not been presented in the literature before.

Materials and Methods
In this research, the research object is the process of propane-butane fraction pyrolysis under the conditions of changing feedstock composition and thermodynamic conditions. The reaction network was developed based on the hydrocarbons reaction network that considers the following reaction occurrence: C3H6→polymers polymers→coke A mathematical model of the pyrolysis process is a system of exponential algebraic equations for calculating the rate constants of the corresponding reactions [17]: where Ki is the rate constant of the i-th reaction, s −1 ; K0i is the preexponential coefficient, s −1 ; P0 is the pressure, atm.; Ei is the activation energy of the i-th reaction, J/mol; and Т is the process temperature, К.
The other 11 equations are differential. They describe the changes in the concentrations of each component of the hydrocarbon mixture during the pyrolysis process. The general view of the equation for each component is given below: where Ci, Ck are the concentrations of i-th and k-th components, mol/L; τ is the contact time, s; and ai,j is the stoichiometric coefficient of the i-th component in the j-th reaction.
The obtained system of differential and exponential equations is solved by the first order Runge-Kutta method.
The coke concentration is calculated as: where G is the feedstock flowrate, kg/s; D is the diameter of the tube, mm; and δ is the coke layer thickness, mm. Due to the fact that the main goal of this work was to use the obtained mathematical model for regulating the operating mode of the furnace and optimizing the process with constantly changing feedstock composition, the process was assumed to be isothermal in order to increase the calculations speed.
The input data used to perform the calculations on the model are presented below. According to the data provided by one of the petrochemical plants in Russia, the pressure drop between the two ends of the pyrolysis coil, with the maximum possible thickness of the coke layer before an emergency stop, is 0.265 MPa.
Thus, we can say that the difference is caused by the pressure loss due to friction against the coil walls and the coke layer along its length: where ΔP is the pressure drop, MPa; λ is the friction loss coefficient along the length; L is the coil length, m; de is the effective diameter of the coil, m; ω is the flow rate of the hydrocarbon mixture, m/s; and ρ is the density of the mixture, kg/m 3 .
According to the technical data provided by the same plant, the heat exchange area is 176 m 2 . However, this is the total value for four parallel coils, so with a coil inner diameter of 100 mm: where Sheat is the heat exchange area, m 2 ; and d is the coil inner diameter, m.
Whereas the cross-stional area of the coil is: The density in this case was calculated using the transformation of the Mendeleev-Clapeyron equation. In this case, the arithmetic mean pressure of the flow at the inlet and outlet of the coil was taken as the pressure: From here the thickness of the coke deposits can be calculated as: where δ is the thickness of coke deposits, mm. Then, the deposition rate during the working cycle between the coke burning stops will be:

Results and Discussion
As part of the work presented in this article, a mathematical model of the reactor was adopted [17]. The adequacy of the adopted model without considering the coke formation was presented by the authors and showed an error of 0.4-0.7 wt%. In the present research, the model was supplemented with equations for formation of polymers and coke, as well as the equations for the rate of coke deposition on the inner surface of the coil. According to the pyrolysis unit at one of the Russian enterprises, the pressure drop in the coil, at which the coke is burned out, is 0.265 MPa. For the given technological parameters, such a pressure drop can be obtained with an average coke layer thickness of about 1 mm (due to an increase in hydraulic resistance in the converging stion of the coil). The average duration of the cycle of the pyrolysis furnace between burns is 1000 h. Thus, the approximate coke deposition rate is 0.0114 mm/day. However, it should be emphasized that this value is an average value over the entire operating cycle.
To verify the model, it is necessary to compare the model-calculated rate of coke deposition in the coil with that calculated from the pressure drop (Equations (4)- (11)) and according to the coke burnout data (Equations (12)- (17)). In this work, it was assumed that the main component of coke is coronene, so that the coke burning corresponds to the following reaction: To burn coke from pyrolysis coil furnaces, the technical air is supplied, and its flow rate varies during the process, as shown in Table 1.
This means, according to the ratio, there is 0.1 kmol in the coke reaction. Accordingly, the mass of coke deposited during the inter-regeneration cycle is: With a coke density of 1493 kg/m 3 , the volume of deposits is: Then, with a pyrolysis coil length of 140 m and diameter of 100 mm, the thickness of coke deposits accumulated during the inter-regeneration cycle is:  Figure 1 compares the model-calculated rate of coke deposition in the coil with that calculated from the pressure drop (Equations (4)- (11)) and according on the coke burnout data.  Figure 1 shows graphs of the coke layer thickness for one day along the length of the pyrolysis furnace coil, expressed by the contact time. The data on the coke layer thickness were calculated by three methods: using the considered model, using calculation methods through an assessment of the pressure loss, and also using the experimental data on the technical air supplied for coke burning. In this case, it is impossible to establish the distribution of thickness along the length of the coil in the last two ways, therefore, they are presented in the form of straight lines.
In general, from the results obtained, we can conclude that the model describes this process quite accurately for the first approximation.
The other way to assess the mathematical model reliability is to compare the calculated amount of coke accumulated in the coil during the whole cycle (42 days) with that calculated using the industrial data on coke burnout. The weight of coke determined using the experimental data was 30.32 kg (Equation (14)). The model calculation results are presented in Figure 2.
This is in a good agreement with value of 30.32 kg calculated using the experimental data. The developed pyrolysis model was tested using the set of input data presented in Table 2. The calculation results are presented in Figure 3. Thus, it can be seen that there is an increase in the concentrations of ethylene and propylene initially, since the high-rate target reactions prevail. This is explained by the fact that the rates of side reactions depend on the concentration of the target products and increase as the target products accumulate. Consequently, side reactions, the formation of polymers and coke, begin to increase towards the end of the length of the pyrolysis furnace coil. As it can be seen in the list of reactions above, the reactants for these reactions are ethylene and propylene, i.e., they are intermediates.
Therefore, approximately 0.5 s after the start of the process, a decrease in the concentration of the target components is observed.
Despite the astronomical time course incorporated in the model, no changes will occur in the dynamics of change in the hydrocarbons' concentration at any time. This is because the only component of this system that changes its value is the thickness of the coke layer. From this point of view, at this stage of the work, it would be more correct to call the model pseudo-steady. Nevertheless, it should be noted that in future work there is a plan to create a mechanism for changing the concentrations of all the components in real time, as well as to create a functional for predicting the state of the system, based on the data accumulated during the installation operation. Figure 4 shows the dynamics of the coke layer growth inside the coil of the pyrolysis furnace during a short period. As presented in Figure 4, the thickness of the coke deposits grows evenly at the beginning of the coil and increases sharply towards the end of the coil through a contact time of 0.4-0.5 s. This happens for the reasons described above, i.e., due to the intensification of side reactions yielding coke. The coke deposition rate in this example is quite high, 0.01 mm per day at its highest point. However, the growth profile is the same every day, whereas in reality the distribution of coke thickness along the pipe is more uneven. The model makes it possible to consider changes in the raw materials composition and technological modes when calculating the rate of formation and deposition of coke.
A slightly different situation is observed if the contact time is reduced ( Figure 5). As can be seen from this graph, the concentrations of ethylene and propylene in this case increase during the entire process time, reaching a certain asymptote, in contrast to the previous case, when both values, reaching some maximum, decreased by the end of the process. Thus, one can suggest that within a reduced contact time period, neither ethylene nor propylene is able to form a significant amount of coke to deposit. Figure 6 shows the dynamics of coke deposition under such conditions. As it can be seen, the absolute value of the coke layer thickness is obviously less, being only about 0.003 mm per day. On the other hand, its growth rate slightly increases: The thickness of deposition begins to grow rapidly 0.2 s into process. Thus, the steadiness in the concentrations of ethylene and propylene is explained, and some amount of these alkenes converts into the coke, but the amount of formed alkenes and that converted to coke per time unit are approximately equal. However, the distribution of deposits along the length of the coil remains the same. This is due to the very same reasons mentioned above in the paragraph concerning coke deposition with longer contact time. If the contact time is increased, the following situation is observed (Figure 7). In this case, the concentration peak of the target components is shifted to the beginning of the coil, while the output concentration is the lowest of the three considered cases. Thus, it makes these conditions the worst of the three. It can be assumed that the hydrocarbon mixture stays too long in the reactor, so side processes, the very same as described for the simulation with 1.2 s contact time, increase their rate and intensity so much so that the output concentrations of target alkenes are lower than in the feed.
The growth of the coke layer thickness, in this case, shows opposite dynamics to that of the previous one. That is, with an increased contact time, the absolute value of the coke layer thickness is the largest in comparison with that of the previous examples, being about 0.07 mm per day. However, the point of coke deposition rapid growth shifted slightly towards the end, at 0.6 s into the process. (Figure 8). In addition, a study was also performed on the effect of the propane/butane ratio in the feedstock on the yield of propylene and ethylene. The results are presented in Figure 9. With an increase in the propane/butane ratio from 4/1 to 1/4 mol/mol, the ethylene concentration decreased from 3.45 mol/L to 3.35 mol/L. Therefore, it is necessary to regulate the composition of the supplied feedstock to obtain the appropriate products in the required quantities.

Conclusions
The presented research yielded a number of important results: (1) The unsteady-state mathematical modeling of hydrocarbon feedstock pyrolysis was developed.
It considers the process of coke deposition along the length of the coil showing that the thickness of the coke layer increases towards the end of the coil. This proves the uneven distribution of coke inside the coil, which is consistent with theoretical concepts of coke deposition processes. (2) The rates of coke deposition along the length of the coil were determined, considering the technological parameters and the composition of the supplied raw materials (the calculated value of coke deposition rate equaled 0.01 mm/day). (3) The influence of the propane/butane ratio in the feedstock on the quality of the products was determined. With an increase in the propane/butane ratio from 4/1 to 1/4 mol/mol, the ethylene concentration decreased from 3.45 mol/L to 3.35 mol/L.
In future work there is a plan to increase the number of components involved in the reaction network to develop mechanisms that simulate the dynamics of the components in real time. This will make the model completely unsteady. There is also a plan to introduce the heat transfer equation into the general system of differential equations, taking into account the rate of heat transfer through the contaminated wall, thereby considering the uneven heating of the flow. In addition, there is as plan to add equations for changing the pressure of the gas flow along the length of the coil by increasing the hydraulic resistance. A model that predicts the state of the system and performs calculations based on the data accumulated for the entire production time is going to be developed. Funding: This research was funded by Russian State Project "Science" FSWW-2020-0011.