Nonsteady-state mathematical modelling of H2SO4-catalysed alkylation of isobutane with alkenes

H2SO4-catalysed isobutane alkylation with alkenes is an important industrial process used to obtain high-octane alkylate. In this process, the concentration of H2SO4 is one of the main parameters. For alkylation, sulphuric acid containing 88%–98% monohydrate is typically used. However, only a H2SO4 concentration of 95%–96% enables alkylate with the maximum octane number to be obtained. Changes in H2SO4 concentration due to decontamination are the main cause of process variations. Therefore, it is necessary to maintain the reactor acid concentration at a constant level by regulating the supply of fresh catalyst and pumping out any spent acid. The main reasons for the decrease in the H2SO4 concentration are accumulation of high-molecular organic compounds and dilution by water. One way to improve and predict unsteady alkylation processes is to develop a mathematical model that considers catalyst deactivation. In the present work, the formation reactions of undesired substances were used in the description of the alkylation process, indicating the sensitivity of the prediction to H2SO4 activity variations. This was used for calculation the optimal technological modes ensuring the maximum selectivity and stability of the chemical–technological system under varying hydrocarbon feedstock compositions.


Introduction
In context of tightening motor fuel requirements, particularly in relation to aromatic hydrocarbon content, there is an increasing demand for alkylate, i.e., the product of the alkylation process of isobutane with alkenes. The obtained alkylate consists of branched hydrocarbons that have a high octane number [1][2][3][4].
In industry, the alkylation of isobutane with alkenes is mainly performed using two types of catalyst, i.e., sulphuric or hydrofluoric acids [5][6][7]. Currently, developments are being made in the field of solid phase alkylation [8][9][10][11][12][13]; nevertheless, to date, only 1% of alkylate is produced using solid catalysts. In 1996, most industrial alkylation plants used Hydrogen Fluoride (HF) as a catalyst; however, owing to its high toxicity, technologies that used a relatively safe analogue, i.e., sulphuric acid, were subsequently developed. Anhydrous HF is a corrosive and highly toxic liquid with a boiling point close to room temperature. Therefore, refineries with HF alkylation plants are under pressure to install expensive mitigation systems minimizing the dangers of HF leaks. Moreover, authorities in many industrialized countries have ceased to license new HF alkylation plants. The other reason for using sulphuric acid is that alkylation gasoline may be contaminated by HF and fluorinated organic compounds that also are highly toxic. Despite other concerns such as the high operating costs, sulphuric acidcatalyzed processes are less hazardous. Sulfuric acid is also a corrosive liquid, but not so volatile, making its handling easier [14].
Unlike alkylation on solid catalysts, the main reason of liquid acidic catalysts deactivation is formation of Acid Soluble Oils (ASO) [2]. The ASO formation mechanism includes olefin oligomerization, disproportionation and cyclization reactions in acid media [10,15]. Accumulation of too much water causes the deactivation of the moisture-sensitive ionic liquid [14]. Effects of aromatics on ionic liquids for C4 alkylation reaction were discussed in [16].
Much work has been devoted to modelling of the solidcatalyst deactivation process [17][18][19][20][21][22][23]. This process and the problems arising from it are influenced by many factors that can be most effectively solved by mathematical modelling [24][25][26][27][28]. It is impossible to completely prevent deactivation; therefore, the only possible way to solve this problem is to create an effective technology for catalyst operation both in the feed cycle (during operation) and in the regeneration process.
Currently, mainly empirical models are used to describe the dynamics of catalyst deactivation. They are not suitable for solving the most important problem of the process efficiency decrease caused by catalyst deactivation. This problem is associated with 90% of all costs for the development and operation of the industrial processes relating to oil refining and petroleum chemistry. Prediction of the catalyst deactivation rate depending on its operating mode is possible only based on a detailed study of the kinetic laws of the deactivation dynamics during the catalytic process. Extending catalyst life by optimising operating conditions and regeneration is extremely important. Revealing the kinetic regularities of the dynamics of the deactivation of liquidphase catalysts, such as sulphuric acid, in the alkylation process makes it possible to predict the occurrence of emergencies and make the necessary decisions for their prevention.
Nevertheless, the methods used for predicting the dynamics of catalyst deactivation are based on empirical approaches rather than on detailed mechanisms and kinetics, which does not allow for the prediction of catalytic system behaviour under changing technological parameters and feedstock compositions.
Several works describe mathematical models of liquid alkylation catalysts deactivation. Modeling of side reactions of isobutane alkylation with butenes catalyzed by trifluoromethane sulfonic acid is described in [10]. Rate equations are first order in catalyst concentration and its Hammett proton activity controlled by ASO content. In [29], the formation pathways for the polymerization of C4 alkenes were discussed and the isomerization between butenes was fully considered. Nevertheless, the proposed kinetic model considers only the detailed mechanism of side reaction, leaving dependence of the reaction rates on the ASO concentration undiscussed.
To date, there is no unsteady mathematical model of H 2 SO 4 -catalysed alkylation of isobutane with alkenes based on industrial data [30].
Models of benzene alkylation with ethylene [31,32] and benzene alkylation with higher alkenes [33,34] have been developed. These models consider the processes of instability and catalyst deactivation by heavy aromatic compounds. For the alkylation of isobutane with alkenes, a mathematical model has been developed that considers the reactions of trimethylpentanes, dimethylhexanes, 2methyldecane, 2,7-dimethyloctane and 2-methylheptane formation [35]; however, it is important to also consider deactivation of the H 2 SO 4 catalyst by acid soluble oils.
The aim of this study was to improve the mathematical model of the H 2 SO 4 -catalysed alkylation of isobutane with alkenes and to describe the mechanism of catalyst deactivation by high molecular weight organic compounds of a polymeric nature in order to consider the emerging deactivation of the catalyst.
In future, the established regularities will be considered in the equations of formal kinetics that describe the rates of chemical reactions occurring in a catalytic reactor [36,37].

Technological scheme of the H 2 SO 4 -catalysed alkylation of isobutane with alkenes
The main components of the system for H 2 SO 4 -catalysed alkylation of isobutane with alkenes are the reactor unit and the crude alkylate separation unit. A schematic block diagram of the process is shown in Figure 1. The plant feedstocks were the Butane-Butylene Fraction (BBF) resulting from the catalytic cracking units and the isobutane fraction from the gas fractionation unit. The mixed raw materials were fed into the contactors of the reactor block via four parallel streams. Each of the four contactors was a horizontal pressure vessel. The piping of the contactors was parallel for hydrocarbons and sequential for the acid. The crude alkylate separation unit included depropanisation, deisobutanisation and debutanisation columns.

System monitoring of H 2 SO 4 -catalysed alkylation of isobutane with alkenes
The main control parameters of the alkylation of isobutane with alkenes were inlet temperature of the feedstock, supply of isobutane coolant, flow rate and concentration of H 2 SO 4 catalyst, feedstock flow rate and circulating isobutane flow rate. The compositions of the BBF, circulating isobutane and isobutane coolant were obtained using gas chromatography during a pilot run of the system for H 2 SO 4 -catalysed alkylation of isobutane with alkenes.
The hydrocarbon composition of the feedstock for the alkylation process is determined by gas chromatography in accordance with GOST 10679-2019. This standard applies to liquefied petroleum gases (propane, propene, butanes, butenes with admixtures of accompanying components and/or their mixtures) obtained during processing of oil, gas condensate, associated petroleum and natural gases. For the analysis, a laboratory gas chromatograph was used, equipped with: (a) a thermal conductivity detector or a flame ionization detector, providing a detection limit for the mass or molar fraction of components of 0.001% or less; (b) a packed or capillary column providing satisfactory separation of the components of the analyzed sample under the conditions of measurements; (c) a column thermostat that provides programmable regulation of the rate of temperature rise and/or maintenance of the set temperature with the limit of the permissible deviation of the thermostat temperature from the average value of no more than 0.1°C in the operating temperature range. Temperature control during measurements should ensure the repeatability of the retention time of the components of one sample no more than 3 s; (d) software that performs the functions of controlling the chromatograph and processing chromatographic data, providing for automated and manual processing of chromatograms, which allows obtaining a graphical image of the chromatogram, identifying peaks, determining the area of the chromatographic peak, calculating the results of analysis by the method of internal normalization or absolute calibration with subsequent normalization of the obtained values, as well as the presentation and storage of data; (e) electronic means of maintaining the speed and pressure of the carrier gas, hydrogen and air flows, ensuring stable retention characteristics of the determined components; (f) a packed or capillary evaporator; (g) a sample introduction system that provides a constant sample volume input from 0.05 to 1.0 mm, depending on the selected instrument configuration. A gas sample dosing system, including a dosing valve made of a corrosion-resistant material, must ensure a repeatability of the injected sample amount of no more than 2% at the required sensitivity to detect the minimum content of the defined component. In this case, the content of components in the measured sample should not exceed the upper linearity limit of the detector; (h) any type of thermostat cooling system (if necessary); (i) valveswitch for backflush (if necessary).
The method for measuring the component composition is based on the separation of sample components using gas chromatography and their subsequent registration with a thermal conductivity detector or a flame ionization detector. The values of the mass or molar fraction of components are determined by processing the chromatographic data by the method of internal normalization using relative mass or molar correction factors or by absolute calibration using absolute molar calibration factors and subsequent normalization of the measured values.
The results are presented in Table 1.
From the point of the product quality, the most favourable components of the process are the butylenes. Therefore, the content of butylenes in the BBF should be maintained at 40-50 vol.% [1]. Simultaneously, the content of diene hydrocarbons should be minimal as these components lead to the formation of undesired complex products when interacting with the acid, thereby diluting it and reducing alkylate quality.
The weight fractions of the alkylate components were determined in accordance with ASTM D5134 using capillary gas chromatography. Equipment used: gas chromatograph with programming the temperature of the column thermostat from 35°C to 200°C, with a heating rate of 1°C/min; capillary column 50 m long, 0.2 mm inner diameter with a grafted (cross-linked) methyl silicone phase, the film thickness of which is 0.5 lm. The procedure for transferring a representative sample of hydrocarbon fluid from the sampling source to a floating piston cylinder was carried out in accordance with ASTM D3700.
The alkylate compositions obtained during an industrial run of the alkylation system are presented in Table 2.
According to the obtained data, during the pilot industrial run of the system, the alkylate contained about 71-74 wt.% of C 8 isoalkanes, providing high octane numbers. About 30-33 wt.% of the C 8 isoalkanes was 2,2,4trimethylpentane.
One of the most important controlling parameters is the isobutane:alkenes ratio. With an excess of isobutane in the reaction zone, target reactions are intensified and side reactions of the process are suppressed. In the current work, the isobutane:alkenes ratio was 11-12 (mol/mol). A significant increase in this parameter is impractical owing to the consequent increase in capital and operating costs.
Considering the large number of parameters affecting the performance of the catalytic alkylation processes, it was necessary to use an integrated approach to mathematical modelling that includes both laboratory and numerical studies. This ensures the consideration of (i) feedstock characteristics and composition, (ii) thermodynamic and kinetic reaction laws in the catalytic reactors and (iii) patterns in catalyst activity change. The mathematical description stems from the sequential performance of the following stages: development of the (i) reaction network, (ii) kinetic equations system, (iii) steady-state reactor model and (iv) unsteady model considering the total volume of raw materials processed and the catalyst activity.

Developing the alkylation process network
On the basis of the performed thermodynamic analysis and experimental data analysis, a formalized reaction network of the H 2 SO 4 -catalysed alkylation of isobutane with alkenes was proposed (Fig. 2). This reaction network included the stages of carbocation formation, reactions of primary alkylation, hydrogen transfer reactions ("self-alkylation"), destructive alkylation, reactions of b-decay, reactions of alkene polymerisation and contamination of H 2 SO 4 with high molecular weight compounds and water.
Based on the proposed reaction network, a kinetic model of the process of H 2 SO 4 -catalysed isobutane alkylation was developed. Expressions for the chemical reactions rates were written in accordance with the law of effective masses.

-catalysed alkylation of isobutane with alkenes
Deactivation of sulphuric acid owing to its dilution with high-molecular compounds and water leads to a decrease in the concentration of protons in the system, which in turn leads to a decrease in the concentration of carbocations (strong Lewis acids), which are intermediate compounds in the target isooctanes formation. The formation of high molecular weight organic substances is due to the cationic polymerisation of alkenes in a H 2 SO 4 medium. In this case, initiation is carried out via the formation of a carbocation from a monomer (alkene) and H 2 SO 4 : The chain growth is an electrophilic addition reaction: The chain growth is an electrophilic reaction. Chain termination can occur due to the addition of an anion of another nucleophile: Chain termination can also be due to proton separation: These reactions are components of the kinetic scheme, along with the formation reactions of high molecular weight hydrocarbons and water, and provide the model's sensitivity to H 2 SO 4 deactivation.
In addition to entering the reaction zone with the feedstock, water in the isobutane alkylation process can also be formed in the following reaction: The spent sulphuric acid can contain dialkylsulfates, some of which can decompose to form oxygen: Via interaction with hydrocarbons, oxygen also forms water. The dilution of H 2 SO 4 with water significantly reduces its catalytic activity; therefore, it is important to maintain the optimal technological mode that ensures minimal water formation. In addition, drying of the butanebutylene and isobutane fractions is required.
The system of heat and mass transfer equations contains a parameter characterising the change in catalyst activity a j. This parameter is determined according to the by-product content and depends on the acid flow rate for regeneration: Initial and boundary conditions: Where C i concentration of the ith component, mol/m 3 ; i = 1,..., N; j = 1,..., M; N: number of substances; M: number of reactions; W j : rate of jth reaction, mol/(m 3 s); s: contact time, sec.; T: temperature, K; q: density, mol/m 3 ; Q j : heat effect of jth reaction, J/mol; C mix p : heat capacity of the reaction mixture, J/(mol Á K); a j : catalyst activity, rel. units.
where k j,initial : initial reaction rate of the jth reaction, k j,current : current reaction rate of the jth reaction.
For the program code of the mathematical model, the Pascal programming language was used in the Borland Delphi 7 programming environment. The choice of this programming language was dictated by the simplicity of the syntax and basic concepts, strict typing, academic character of the language and the availability of a license for the corresponding software.
To simulate the alkylation process, the following data were used: 1. composition of raw materials and products and their characteristics (density, heat capacity, viscosity, etc.); 2. technological parameters (temperature, pressure, flow rates); 3. characteristics of the catalyst; 4. equipment characteristics (geometrical dimensions, configuration, dimensions of branch pipes, etc.); 5. technological regulations and technological scheme.
The data for the process modeling were taken from the factory databases, which were obtained on the existing equipment in the factory using a telemetry system.
The inverse kinetic problem was solved by the following way. To calculate the activation energy and the pre-exponential factor in the Arrhenius equation, we applied the theory of absolute velocities of chemical reactions (the theory of a transition state).
The pre-exponential factor in the Arrhenius equation was calculated using the formula: where v: transmission coefficient; k: Boltzmann's constant (1.38 Â 10 À23 J/K); T: temperature, K; h: Planck constant (6.62 Â 10 À34 J/s), DS # : activation entropy, kJ/mole, R: universal gas constant (8.314 kJ/(mole K). The activation energy was calculated as where DH # : activation enthalpy, kJ/mole. The search for the transition state was performed by the method of quadratic synchronous transit QST2.

Verification of the developed mathematical model
The developed model of H 2 SO 4 -catalysed isobutane alkylation was software implemented in a high-level language. The developed model of H 2 SO 4 -catalysed isobutane alkylation was software implemented in a high-level language in the Borland Delphi 7 programming environment. By solving the inverse kinetic problem and based on the data obtained during the experimental industrial run of the system, the kinetic parameters of the chemical reactions were determined and used in a kinetic description of the process. The model was verified by comparing the calculated and experimental alkylate compositions (Fig. 3).
According to the results presented in Figure 3, the model showed satisfactory convergence when calculating the commercial alkylate composition. The average calculation error was less than 2 wt.%.

Assessment of the main control parameters influencing alkylate yield and quality
With the developed model, an assessment of the main control parameters influencing the alkylate yield and quality was performed (Figs. 4 and 5). The octane number of the mixture can be represented as the sum of two components: additive and non-additive: where ON: octane number of gasoline mixture; B: total deviation of octane numbers from additivity; C i : concentration of the ith component, rel. units.
where B i , B j : units characterizing the tendency of the ith molecule to intermolecular interaction with the jth molecule, which can be expressed through the dipole moments of the molecules: where a and g: kinetic parameters that determine the intensity of intermolecular interactions depending on the dipole moment D, for a toluene molecule are numerically equal to 2.21 and 1.09, respectively; D max : maximum dipole moment of C 9+ aromatic hydrocarbon molecules [38]. Changes in the load and composition of the raw material, the activity of the catalyst (sulfuric acid), and the technological mode parameters give rise to an area of uncertainty in determining the octane number and reaction rate constants, which are known within certain intervals. The size of the uncertainty region (for an octane number of 1-1.5 units, Fig. 4) is additionally affected by the range of quantity, type, and availability of the source data. The uncertainty region determines the range of parameter changes and the margin factor for the optimized characteristics, including the octane number. Increase of the feedstock initial temperature increased the rates of both the target and side reactions; therefore, there was an optimal temperature range that depended on many factors such as composition of the feedstock, acid concentration, isobutane:alkenes ratio and configuration of the technological equipment.
Increase of the circulating isobutane flow rate increased the isobutane:alkenes ratio, which favourably affected alkylate quality and yield (Fig. 6); however, it is impractical to maintain significantly high values of this parameter owing to the consequential high capital and operating costs.
An increase in the supply of circulating isobutane first leads to an increase in the yield of alkylbenzene due to an increase in the degree of conversion of the initial butylenes, as well as suppression of side reactions of the polymers formation. Then, the yield curve of alkylbenzene reaches steady state.
The optimal consumption of circulating isobutane is determined based on the content of butylenes and isobutane in the butane-butylene fraction stream (process feedstock) and in the circulating isobutane stream itself. The calculation results for two different compositions of the butane-butylene fraction (all other things being equal for the process of sulfuric acid alkylation) are shown in Figure 7.
The total content of butylenes is 47.9 and 44.8 wt.% for compositions of raw materials No. 1 and No. 2, respectively. Thus, in order to achieve the maximum yield of alkylate, when processing raw materials close in composition to feedstock No. 1, the optimal circulating isobutane flow rate lies in the range 45-50 m 3 /h, and when processing raw materials close in composition to feedstock No. 2, the optimal circulating isobutane flow rate lies in the range of 37-42 m 3 /h.
A study of the influence of the butylene content in the butane-butylene fraction on alkylate yield was carried out under constant technological parameters (Tab. 3).
The calculation results are presented in Figure 8. According to the results presented in Figure 8, an increase in the BBF butylenes content led to increases in the alkylate yield and its octane number. This was due to intensification of the target reactions of the formation of trimethyl-substituted octane isomers.
Calculations were also carried out on a model for various operating modes of the isobutane alkylation process. The technological parameters are presented in Table 4. H 2 SO 4 concentration was one of the main process parameters. For alkylation, sulphuric acid containing 88-98% monohydrate is typically used. However, only 95-96% H 2 SO 4 allows for alkylate with the maximum octane number to be obtained [35]. Changes in H 2 SO 4 concentration due to decontamination are the main reason for process instability; therefore, it is necessary to maintain the acid concentration in the reactor at a constant level by regulating the supply of fresh catalyst and pumping out the   spent acid. The main reasons for the decrease in H 2 SO 4 concentration are the accumulation of high molecular weight organic compounds and dilution by water. The occurrence of irreversible side reactions under the thermobaric conditions of the experiment was confirmed by the thermodynamic analysis performed. The kinetic regularities of these reactions were established by solving the inverse kinetic problem using the experimental data.
The mathematical model of the alkylation reactor includes the functional dependence of the liquid-phase catalyst activity change on the amount of by-products formed, allowing us to calculate the fresh acid flow rate required to ensure maximum selectivity and stability of the process conditions under conditions of changing feedstock composition. According to data presented in Table 4, the parameters of the H 2 SO 4 -catalysed alkylation system for isobutane varied over wide ranges.
The catalyst activity in the alkylation process was determined using the concentration of H 2 SO 4 fed to the contactors. To maintain the concentration of waste acid at the level of 90 wt.%, fresh acid was added to the system as an appropriate volume of spent acid was removed. A graph of the catalyst change activity during the alkylation process is shown in Figure 9.
According to the data presented in Figure 7, the catalyst activity at the outlet of the reactor block was 94.3 and 93.8 wt.% for modes 5 and 7, respectively; whereas for the other modes, the catalyst activity was 95-96 wt.%. An increase in the side reactions rate, in particular the formation of high molecular weight hydrocarbons during cationic polymerisation reactions in modes 5 and 7, can be explained by the relatively low isobutane:alkenes ratios of 9.2 and 9.4 (mol/mol), respectively. For the other modes, this ratio was maintained at a level of 11-12 (mol/mol), which made it possible to suppress undesirable alkene polymerisation reactions.
The mathematical model included the functional dependence of changes in the liquid-phase catalyst activity on the amount of by-products formed. This allowed calculation of the fresh acid flow rate required to ensure maximum selectivity and stability of the chemical-technological system under conditions of changing hydrocarbon feedstock composition.
The optimal consumption of fresh acid is determined based on the condition that the concentration of the spent acid for the first stage of the contactors should be at least 93-95 wt.% since then this acid sequentially enters the second stage, and if its concentration drops below 90 wt.%, intensification of side reactions of alkenes polymerization will occur. With an increase in the BBF supply, there is an increase in the acid consumption and, as a consequence, a decrease in its concentration in the contactors. The calculation results are presented in Figure 10.
With an increase in the BBF flowrate, an increase in the consumption of sulfuric acid occurs. Therefore, fresh acid feed is required to maintain the acid concentration at the outlet of the contactor at a level of 93-95 wt.%.
In this study, the experimental data reflected changes in the H 2 SO 4 catalyst concentration due to the formation of  high molecular weight hydrocarbons in the isobutane alkylation process. The developed mathematical model of the alkylation process was sensitive to changes in feedstock composition and showed a high predictive potential in relation to alkylate yield and quality. An example of the optimum process conditions for a sulfuric acid alkylation is presented below (Tabs. 5 and 6). These technological parameters provide a high yield of alkylate (57.5 kmol/h per one contactor) and a high octane number of alkylate (95 points according to the research method).
However, these recommendations are suitable only if the processed feedstock composition is close to the presented above. When processing feedstock of a different composition, the optimal parameters will be different. In addition, when changing the loading of the plant for BBF, the optimal flow rate of circulating isobutane and the temperature in the contactors would also be adjusted.

Conclusion
The current research is devoted to developing the unsteady mathematical model of H 2 SO 4 -catalysed alkylation of isobutane with alkenes and to optimizing the technological modes. The study yielded a number of findings.
The proposed level of reaction network detalization considers the high-molecular compounds formation and the process unsteadiness. This ensures the universality and adequacy of the reactions kinetic description and preserves the mathematical model sensitivity to the hydrocarbon feedstock composition. Considering reactions of high molecular weight hydrocarbons and water formation provides the model's sensitivity to H 2 SO 4 deactivation and allows revealing the optimal technological modes. Thus, for the feedstock containing 0.3 wt.% of butadiene-1,3, 15.75 wt.% of butylene, 43.41 wt.% of isobutane, 5.76 wt.% of iso-butylene, 8.43 wt.% of n-butane, 16.07 wt.% of cis-butylene, 9.92 wt.% of trans-butylene, the optimal H 2 SO 4 flowrate to contactors 1 and 2 is 3.20 m 3 /h and 3.58 m 3 /h, respectively.
Further research will aim to complement the developed kinetic model with heat balance equations and calibrate the kinetic parameters of the model over a wider range of experimental data.