## 1. Introduction

Hydrocarbon steam cracking is the most important process for producing industrial chemicals such as olefin and aromatics. Mathematical models for steam cracking simulation have been studied for several decades and various models have been developed such as SPYRO [1] and COILSIM [2]. Meanwhile, the steam cracking model can be applied to the steam furnace design and optimization [3]. As various feedstock and furnaces are used in steam cracking, a more accurate and robust model is needed.

A mathematical model is a collection of relations among variables representing certain properties of a system, using equations and inequalities. By solving a model, one can predict the values of some properties of the system, given the values of some other properties, or predict the distributions of the values of some properties in particular domains, given certain boundary conditions. To meet the scientific and engineering demands in terms of revealing the characteristics of a system in-depth, a mathematical model may involve relations of properties at different scales of the system [4], referred to as a multi-scale model. The multi-scale model often consists of nonlinear equations and differential equations, and is not easy to solve. The information communication between scales is the key factor in multi-scale models [5]. In recent years, multi-scale modelling has been applied in all fields of chemical engineering such as thermodynamics, reaction engineering, polymer materials and CFD (computational fluid dynamics), among others. Multi-scale modelling is used in steam cracking models to reveal the nature of steam cracking and to establish a more accurate and robust model.

The kernel part of the steam cracking model is the reaction network. Researchers have developed various reaction models to describe the steam cracking process. There are three different types of reaction models: empirical models, molecular models and elementary reaction models. Empirical models use a large database of experimental results to regress a number of empirical correlations for the yields of the main products as a function of a number of easily measurable process variables [6]. Empirical reaction models need a large number of experiment data to render a regression. Once the feedstock or the furnaces change, the reaction model has to be re-established to obtain accurate simulation results. On the basis of empirical reaction models, molecular reaction models have been developed and frequently used [7-9], e.g., the Kumar reaction model [10]; 22 molecular reactions are involved in Kumar’s reaction model: one primary reaction and 21 secondary reactions. The primary reaction is shown as eq. (1). The selective coefficients of the primary reaction (*a*_{1}*~a*_{10}) have been regressed from experimental data. If the feedstock or operating conditions change, the selective coefficients of the primary reaction should be regressed again.

Nowadays, elementary reaction models have been widely used to develop a more reliable and robust mathematical model. Since the pioneering work of Rice [11-13] there has been a general consensus about the elementary reaction mechanism. As the mathematical difficulties encountered for solving the detailed kinetic models can be overcome by the development of stiff solvers [14], detailed elementary reaction networks are widely used to accurately describe chemical processes in a wide range of process conditions and feedstock. Sundaram et al. [15] established a radical reaction model for the pyrolysis of simple paraffins, olefins and their mixtures. Scharfe et al. [16] established a radical reaction model for n-hexane pyrolysis. Joo et al. [17] established a radical reaction model for industrial naphtha cracking furnaces.

However, it has remained problematic to generate thousands of elementary reactions and to determine the reaction rate constants of these reactions until an automatic reaction network generation technique was studied. Today, computers are used not only to solve the simulation numerically, but also to generate the elementary reaction network, construct the model and calculate the kinetic parameters [2]. Many research groups have developed computer tools for automatically generating these mechanisms [18-23] such as RMG [24] and CRACKSIM [25]. The elementary reaction model is expected to be used in furnace design and operating condition optimization. However, the huge number of reactions in the elementary reaction network usually requires significant CPU time during simulation. Thus, the steam cracking simulation model has to be reduced before it can be used in operating condition optimization, in order to make real-time optimization (RTO) realistic.

There are multiple methods for creating model reduction. In terms of chemistry, quasi-steady-state approximation (QSSA) [26], reaction rate analysis [27] and reaction path analysis [28], among others, can be used to reduce a reaction model. In terms of mathematics, principle component analysis [29] and sensitivity analysis [30] can also be used. On the other hand, data driven methods can also be used in model reduction, for example, the black box model, a neural network [31] and PCA based ROMs [32]. Several assumptions can also be applied to help retain the mechanism at a manageable size. The most important assumption is the *μ* radical hypothesis, which assumes that bimolecular reactions can be neglected for radicals with more than five carbon atoms [33]. The latter are also referred to as *μ* radicals. Thus, the QSSA method can be applied to remove *μ* radicals from the reaction network.

A multi-scale model for the steam cracking process is established in this chapter. The multi-period steam cracking process is also studied in the context of the established multi-scale model. Coking is an unavoidable factor during the multi-period steam cracking process. Coking increases pressure drops in the reaction tube, decreases the coefficient of heat transfer between the furnace and tube, and most importantly, increases the outer-wall temperature of the tubes. If one of the tubes in the furnace reaches the maximum allowance temperature of the tube material, the furnace must be shut down to execute a decoking process, otherwise the tube will be destroyed. The operation conditions are generally maintained constant during the long-term steam cracking process. Thus, dynamic operating conditions need to be optimized using a detailed steam cracking model in order to achieve a higher profit. Abel et al. [34] used the SQP method to solve a real-time optimization problem in the olefins production process. Tarafder et al. [35] proposed a multi-objective optimization problem in the operation and design of a styrene manufacturing process. Li et al. [36] applied an artificial neural network (ANN) hybrid model in the operation optimization of a naphtha industrial cracking furnace. Gao et al. [3] used a new parallel hybrid algorithm combining NSGA-II with SQP on multi-objective optimization for the periodic operation of the naphtha pyrolysis process. However, due to the complexity of the elementary reaction model, the researchers did not use this reaction model in the optimization problem. In this chapter, an elementary reaction model is applied to the operating condition optimization problem to obtain a more reliable optimization result. Based on this, a surrogate coke thickness model is proposed to make multi-period optimization with parallel simulation possible.

The general idea of this chapter is outlined in Figure 1. The first step in the conceptual development of a detailed molecule-based model for a complex feedstock is to determine an accurate molecular representation of the feedstock. Then, a multi-scale steam cracking model is established following the feedstock prediction. Finally, operating condition optimization of multi-period steam cracking is carried out using the established multi-scale model.

This chapter is structured as follows. Section 2 offers a detailed discussion of the establishment of the multi-scale steam cracking model; a case study for naphtha steam cracking simulation is presented. Section 3 provides the operating condition optimization model; surrogate coke thickness model and parallel simulation are used to accelerate the computing of the optimization model. A case study of operating condition optimization is presented and the results are presented and discussed. Finally, section 4 offers conclusions from this study.

## 2. Steam cracking model

### 2.1. Multi-scale model for steam cracking

A multi-scale model is proposed in this chapter to reveal the nature of the steam cracking process and to generate a reliable and robust model for accurately predicting the yields of products. The multi-scale model is built up from the molecular level, to the reaction level to the process level. The establishment of the multi-scale model is discussed in detail in the following sections.

Conventional analytical techniques are generally incapable of directly measuring the identities of all the molecules in complex feedstock, especially for the high carbon number range; however, this applies only to indirect characteristics [37]. Here, the Shannon entropy method [38] and probability density function were used to predict detailed feedstock composition, based on the analytical data. A detailed model for feedstock composition prediction can be found in [39]. The objective function of this NLP (non-linear programming) problem is shown in eq. (2), where *S(x)* is Shannon entropy and *x*_{i} is the mole fraction of component *i*.

Figure 2 shows the predicted feedstock composition using Shannon entropy, compared with the experimental data. The predicted results show the effectiveness of the Shannon entropy theory in obtaining the missing information for models of the steam cracking process.

A process-level model consists of mass balance equations, momentum balance equation and energy balance equation (eqs. (3-5)) [40, 41]. For the simulation of smooth tubular reactor types, the use of a one-dimensional reactor plug-flow model is generally recognized as providing a sufficient degree of accuracy, as all radical profiles are wiped out due to the high turbulence corresponding to Reynolds numbers of over 250 000. The plug-flow reactor model implicitly assumes that there is no mixing in the axial direction, but rather, perfect mixing in the transverse direction [2].

In eq. (3), *N*_{m} is the concentration of species *m* and *L* is the length of the reactor tube. *S* is the flow area, *V* is the volume flow rate, *μ*_{im} is the stoichiometric coefficient of reaction *i* and *r*_{i} is the reaction rate of reaction *i*. In eq. (4), *P* is the pressure of the mass flow, *f* is the Fanning friction factor, *E(L)=L*_{e}*/L*, *L*_{e} is the equivalent length of the reactor tube, *G* is the mass flow rate, *D*_{i} is the internal diameter of the tube and *ρ* is the density of the gas mixture. In eq. (5), *T* is the temperature of the mass flow, *D*_{o} is the outer diameter of the tube, *k* is the overall heat transfer coefficient, *T*_{w} is the outer wall temperature of the tube, *ΔH*_{fm}^{0} is the standard heat of formation of species *m* and *c*_{pm} is the specific heat of species *m*.

It should be noted that the measuring point of the COT (coil-outlet temperature) of an industrial furnace is usually on the outer wall of the tube. Thus, the measured COT has a temperature difference from the outlet temperature of the gas mixture. Eq. (6) has been derived from heat balance equations in order to adjust to the temperature difference; the results show that there exists a 15-20K temperature difference, which agrees well with what has previously been reported.

In eq. (6), *λ*_{isolation} is the heat transfer coefficient of the isolation layer. *T*_{isolation_o} is the outer wall temperature of the isolation layer, *r*_{o} is the external diameter of the tube, *r*_{i} is the interior diameter of the tube, *d*_{isolation} is the thickness of isolation layer, *d*_{coke} is the thickness of coke and *α*_{mixture} is the heat transfer coefficient of the gas mixture in the tube.

The reaction model is the most important part of the steam cracking model. Many researchers have conducted in-depth studies on the elementary reaction model. An elementary reaction model can contain thousands of reactions and hundreds of species. The reaction model can be extremely hard to solve due to its stiffness. An accurate and robust elementary reaction model is developed in this chapter, and a Gear algorithm is used to solve the stiff ODEs. Generally, a detailed reaction network is generated by allowing the feedstock components to react according to different reaction families. The reaction families can be summarized as follows: (1) initiation reaction and termination reaction: *R*_{1}*—R*_{2}
*R*_{1}*∙+R*_{2}*∙*; (2) hydrogen abstraction reaction: *R*_{1}*—H+R*_{2}
*R*_{1}*∙+R*_{2}*—H*; (3) radical addition and *β*-scission reaction: *R*_{1}*∙+R*_{2}*=R*_{3}
*R*_{1}*—R*_{2}*—R*_{3}*∙*.

The hydrocarbon steam cracking elementary reaction network can be separated into several sub-models. based on the composition of feedstock as shown in Figure 3. The sub-models are generated separately, based on the reaction families described above.

An elementary reaction model for light hydrocarbon can be found in much of the literature and databases, and it is more accurate than the automatic generated elementary reaction model. Thus, the elementary model for steam cracking is separated into two parts: a light hydrocarbon and heavy hydrocarbon reaction model (carbon number greater than five). The light hydrocarbon reaction model is generated using RMG [24]. RMG considers each species as unique and comprising a set of molecular structural isomers. When a reaction network is generated using RMG, it need not consider all the isomers in the real steam cracking process; instead, only a set of representative species are considered during the generation of a reaction network. The heavy reaction model is combined using different reaction models of pure compound feedstock and each reaction model is generated from the reaction families. The reaction coefficients can be obtained from the summary of experimental data (Table 1) in the work of Dente et al. [26]. The heavy hydrocarbon reaction model was generated using our own code.

The automatic generated reaction network may contain a large number of unimportant reactions and species. These reactions can increase the complexity of the model and make the model hard to solve. Thus, reaction model reduction is needed following the automatic model generation. The QSSA method is introduced first to remove the *μ* radicals in the reaction network [33]. As eq. (6) shows, the reaction rate of *μ* radicals can be treated as zero, based on the assumption.

Thus, eq. (8) can be derived from eq. (7):

In eq. (8), *r*_{j} is the rate of direct formation of *j*-radical, *k*_{i,j} is the rate constant for the isomerization reaction (*R*_{i}
*R*_{j}) and *k*_{dj} is the total rate constant for the decomposition reactions of *R*_{j}; *μ* radicals are reduced using eq. (8). The number of species included in the model is decreased.

Reaction rate analysis is used to reduce the unimportant reactions in the reaction network. The average reaction rate in the reaction network reflects the importance of the reaction in the network. Thus, we can rank the reactions based on the average reaction rate and reduce the reactions where the reaction rate is less than *R*_{min} [44].

The flowchart of the automatic generation and reduction of the reaction network is shown in Figure 4.

The physical properties of some species (radical, non-common substances, etc.) involved in the model are difficult to obtain from databases. The physical properties of these species can be calculated using the group contribution method [42]. RMG also supplies a thermochemistry estimates utility using the group contribution method and was used in our model to automatically calculate the physical properties of these species.

### 2.2. Case study: Naphtha steam cracking model simulation

A light hydrocarbon reaction model was generated using RMG, which contained 91 reactions and 26 species (Table 2). A naphtha steam cracking reaction network was generated based on the light hydrocarbon network. The naphtha reaction network contained 2424 reactions and 125 species.

Thirteen sets of industrial data (Table 3) were used here to verify the established multi-scale steam cracking model. The industrial data were taken from a steam cracking furnace designed by KBR (Kellogg Brown & Root). The data were collected twice a day and each set of data was the average value of these two parallel experiments as a means for preventing any errors. The multi-scale steam cracking simulation took roughly 70s CPU time. The simulation results and the industrial data of the mass fraction of the main products are shown in Figure 5, where the x-axis represents the industrial data and y-axis represents the simulation results. All points in Figure 5 are around the 1:1 diagonal line, which shows that the error of most results were small. Thus, the established multi-scale model was defined as accurate and robust; it satisfies industrial needs and can be applied further in our study of operating conditions optimization.

## 3. Multi-period steam cracking optimization

### 3.1. Multi-period optimization model for steam cracking

Coking is an unavoidable factor in the long-term steam cracking process. During the steam cracking process, coke forms on the inner walls of the tube. With coke formation, the internal diameter of the tube decreases, pressure drops increase and outer wall temperature increases. It is generally accepted that coke forms from unsaturated hydrocarbons and aromatics [43]; on this basis, an empirical model for coke formation rate is proposed, as shown in eq. (9).

In eq. (9): 1=C_{2}H_{2}; 2=C_{2}H_{4}; 3=C_{3}H_{6}; 4=1-C_{4}H_{8}; 5=C_{4}H_{6}; 6=C_{6}H_{6}; 7=C_{7}H_{8}; 8=xylene; 9=C_{2}H_{3}-C_{6}H_{5}.

Steam cracking is a dynamic process, as coke grows inside the tube. However, coke formation is slow enough that we can divide the entire cracking period into a series of virtual steady state periods. Thus, the established multi-scale model can be used in each steady state period, with coke thickness updated between periods.

Operation condition optimization is carried out based on the multi-period process model. Coil outlet temperature (COT) is selected as the variable to be optimized. Thus, the COT of all time periods of the multi-period model is discretized as [*COT*_{1}*, COT*_{2}*,..., COT*_{n}].

The optimization model is summarized below.

Subject to:

In objective function (10), *M* is the number of considered species, *N* is the period number, *ω*_{t} (*t*=1, 2,..., *M*) are weighted fractions based on the price of each species and *y*^{o}_{j,t} is the mass fraction of selected species *j* in products of period *t*. Eqs. (11-13) describes mass balance, momentum balance and energy balance equations in period *t* (*t*=1, 2,..., *N*). In eq. (14), the internal diameter of period *t*+1 *D*_{i,t+1} equals the internal diameter of the previous period *D*_{i,t} minus the coke layer thickness. In eq. (16), peak outer-wall temperature should not exceed the maximum temperature of the tube material. Eq. (17) shows that the adjacent COT difference is restricted to a certain region to keep the operation stable. Eq. (18) shows the upper and lower boundaries of optimization variables.

The optimization procedure is shown in Figure 6.

### 3.2. Surrogate coke thickness model

The most time-consuming part of the optimization is the simulation of the multi-period model. As Figure 6 shows, the only connection between adjacent periods is the coke thickness. Coke thickness is related to the feedstock component, furnace running time and operating conditions. As it has been assumed that the feedstock component is fixed between periods and only COT changes in operating conditions, as eq. (19) shows, then *d*_{t} is the coke thickness in period *t*.

The accumulated coke thickness within a certain period *δd*_{t} is assumed to only be related to the furnace running time and *COT*_{t} (*t*=1, 2,..., *N*-1), shown as eqs. (20) and (21).

Thus, *δd*_{t} in period *t* can be regressed using the polynomial function shown as eq. (20). Here, the coke thickness data is generated using the original multi-period simulation model and based on this, a surrogate coke thickness model is obtained through regression. The coke thickness using the surrogate model and original multi-period simulation model are shown in Figure 7. Dots in Figure 7 are coke thickness from the original multi-period simulation model and the surface is from surrogate model.

The coke thickness from the surrogate model fits well with the original model; thus, the decoupled multi-period cracking model, combined with the surrogate coke thickness model was used in the multi-period simulation. The initial coke distribution along the serial operation periods was carried out using the surrogate model. Thus, the multi-period simulation problem was decoupled into *N* sub-problems and simulated, respectively, in parallel, as shown in Figure 8.

### 3.3. Optimization with parallel simulation

Figure 8 shows optimization with a parallel simulation procedure using the surrogate coke thickness model.

In Figure 8, a surrogate coke thickness model is regressed from the results of the original multi-period model and the surrogate coke thickness model generates the coke thickness distribution for each period prior to the simulation. Once the surrogate model is regressed, the process for each period can be simulated in parallel. The simulation results are sent to the optimization model. If the criteria are satisfied, the optimization stops; if not, the optimization model generates a new set of COT and returns it to the simulation model.

There are several common types of parallel computing methods: phase parallel, divide and conquer parallel, pipeline parallel, master-slave parallel and work pool parallel methods. In this instance, the work pool parallel job partitioning method was used in the optimization. Simulation for one period can be treated as one job in the parallel simulation. All the jobs are stored in the work pool; processors fetch jobs from the work pool as long as the work pool is not empty. A detailed scheme for parallel simulation is shown in Figure 9.

In Figure 9, *NP* is the resources that can be used in a parallel simulation (*NP*=computer number×processor number for each computer) and *NJ* is the period number in the multi-period simulation. If *NP≥NJ*, that means all jobs can be calculated in one iteration. Otherwise, *NP* jobs must be fetched from the work pool and be postponed until all processors are idle before the next iteration. Simulation will continue until the work pool is empty.

Another benefit of using parallel simulation is a warm start. To accelerate the steam cracking simulation, iteration information is stored after the first simulation, which is called a warm start. If the iteration information is recorded, it can be used as an initial value in the next simulation for reducing CPU time. A comparison between a warm start for serial and parallel simulation is shown in Figure 10.

Iteration information for the most recent period is stored in serial simulation. The iteration information for each period is not entirely the same. Inconsistency of iteration information makes the simulation slow. In a parallel simulation, however, iteration information for each period can be stored. This delivers a faster simulation for each period.

### 3.4. Case study: Naphtha steam cracking model optimization

An industrial multi-period steam cracking instance was used as a case study. Input data is shown in Table 4. The detailed feedstock composition was calculated using the Shannon entropy method [39]. An elementary reaction model and optimization with parallel simulation was used in this case study. The computer used in this optimization had eight cores of i7-2600 CPU and 4GB memory.

The optimization results show a 0.62% ethylene increase compared to the invariant operating condition that was implemented in the practice. The comparison of the serial optimization results and optimization with parallel simulation results are shown in Figure 11 and Figure 12. It can be seen that the tendency of serial and parallel optimization results were the same. There was a high outlet temperature and high conversion of the major products in the beginning and intermediate periods, and the outlet temperature of final periods were rather low for reducing the coke formation and for avoiding a high outer-wall temperature.

The error of parallel simulation compared to serial simulation is shown in Table 5.

The CPU time of serial optimization was 17.78hr, while the CPU time of parallel optimization was 2.08hr. There was a 8.55 x speedup compared with serial optimization, because parallel calculation and a warm start strategy were used. Thus, the operating conditions could be dynamically optimized to track the changing market conditions.

The error in optimization with the parallel simulation model was caused by the surrogate coke thickness model. In eq. (22),

The feedstock composition and operating conditions, except COT, were fixed. However,

Accumulate coke thickness in each period was related to the operating conditions of all previous periods. Eq. (23) is simplified as eq. (22) to make the model easier to regress and the error caused by the simplification is acceptable according to the error shown in Table 5.

## 4. Conclusion

A multi-scale steam cracking model was established in this chapter. An elementary reaction network was generated to obtain a more accurate and robust model. The results showed that the established multi-scale model agrees well with the industrial data. A surrogate coke thickness model is thus proposed to make multi-period optimization with parallel simulation possible. A case study showed that multi-period optimization with parallel computing possesses an 8.55 x speedup compared to optimization with serial simulation. Parallel computing makes real-time optimization (RTO) possible in multi-period optimization. The average error of optimization using a parallel simulation model was 0.04% compared to optimization using a serial model. The error was deemed acceptable for the optimization of large-scale industrial steam cracking processes.

Optimization using the parallel computing method can also be used in other aspects of chemical processes. Due to the complexity of chemical processes, it is always difficult to conduct a simulation or optimization using a multi-scale model. The connections of different levels in a multi-scale model or different equipment in chemical process can be decoupled and transferred to a multiple sub-system model, which can be simulated in parallel. This method can bring powerful computing performance into play in chemical engineering.