Simulation and Optimization of Multi-period Steam Cracking Process

Abstract

Hydrocarbon steam cracking is the most important process for producing industrial chemicals such as olefin and aromatics. Steam cracking modelling and optimization is an effective way for increasing production and saving energy. In this chapter, multi-scale modelling and elementary reaction networks are established and used in the modelling and optimization of steam cracking. However, the large scale of the optimization model makes it difficult to obtain a solution. Thus, a surrogate coke thickness model for long-term steam cracking is proposed in this chapter to remove the connection between different periods of steam cracking process. By so doing, a parallel simulation can be used to accelerate optimization. An industrial case study showed optimization time to be significantly reduced from 17.78 hours to 2.08 hours using multi-period optimization with parallel simulation and the surrogate coke thickness model. It has been shown that a 0.62% increase in ethylene yield can be obtained via operating condition optimization, which demonstrates the effectiveness of the multi-scale steam cracking model and multi-period optimization with parallel simulation.

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₁~a₁₀) have been regressed from experimental data. If the feedstock or operating conditions change, the selective coefficients of the primary reaction should be regressed again.

$$\begin{array}{l}\text{naphtha}\to a_1{\text{H}}_{\text{2}}+a_2{\text{CH}}_{\text{4}}+a_3{\text{C}}_{\text{2}}{\text{H}}_{\text{4}}+a_4{\text{C}}_{\text{2}}{\text{H}}_{\text{6}}+a_5{\text{C}}_{\text{3}}{\text{H}}_{\text{6}}+\\ +a_6{\text{C}}_{\text{3}}{\text{H}}_{\text{8}}+a_7{\text{C}}_{\text{4}}{\text{H}}_{\text{10}}+a_8{\text{C}}_{\text{4}}{\text{H}}_{\text{8}}+a_9{\text{C}}_{\text{4}}{\text{H}}_{\text{6}}+a_{10}{\text{C}}_{\text{4}}\text{'s}\end{array}$$E1

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.

Figure 1.

Diagram of simulation and optimization of the steam cracking process methodology

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.

1. Feedstock composition prediction

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.

$$\max S(x)=-{\displaystyle \sum }_{i=1}^n{x}_i\ln x_i$$E2

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.

Figure 2.

Comparison of experimental and calculated feedstock component mole fractions [39].

1. Process-level model

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].

$$\frac{d{N}_m}{dL}=\frac{S}{V}{\displaystyle \sum _i{\nu }_{i,m}{r}_i}=f_N\left(T,P,N_m\right)$$E3

$$\frac{dP}{dL}=-\frac{f\cdot E(L)\cdot G^2}{5.07\times \rho \cdot D_i\times {10}^4}=f_P\left(T,P,N_m\right)$$E4

$$\frac{dT}{dL}=\frac{k\pi D_o(T_w-T)-{\displaystyle \sum _m\Delta H_{fm}^0\cdot \frac{d{N}_m}{dL}}}{{\displaystyle \sum _m{c}_{pm}{N}_m}}=f_T\left(T,P,N_m\right)$$E5

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⁰ 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.

$$\Delta COT=\frac{{\lambda }_{isolation}\left(COT-T_{isolation\_o}\right)}{\ln \frac{r_o+d_{isolation}}{r_o}}\left(\ln \frac{r_i}{r_i-d_{coke}}+\frac{\ln \frac{r_o}{r_i}}{{\lambda }_{tube}}+\frac{1}{{\alpha }_{mixture}(r_i-d_{coke})}\right)$$E6

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.

1. Reaction-level model

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₁—R₂ $\leftrightarrow$ R₁∙+R₂∙; (2) hydrogen abstraction reaction: R₁—H+R₂ $\leftrightarrow$ R₁∙+R₂—H; (3) radical addition and β-scission reaction: R₁∙+R₂=R₃ $\leftrightarrow$ R₁—R₂—R₃∙.

Figure 3.

Sub-models in the generation of hydrocarbon steam cracking reaction networks.

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.

Initiation reactions: unimolecular decomposition of C-C bonds
CH3-Csec	Csec-Csec	Csec-Cter	Csec-Cquat
5×10^16exp(-83,500/RT)	5×10^16exp(-81,000/RT)	5×10^16exp(-80,000/RT)	5×10^16exp(-78,000/RT)
H-abstraction reactions of alkyl radicals
	Primary H-atom	Secondary H-atom	Tertiary H-atom
Primary radical	10^8exp(-13,500/RT)	10^8exp(-11,200/RT)	10^8exp(-9,000/RT)
Secondary radical	10^8exp(-14,500/RT)	10^8exp(-12,200/RT)	10^8exp(-10,000/RT)
Tertiary radical	10^8exp(-15,000/RT)	10^8exp(-12,700/RT)	10^8exp(-10,500/RT)
Isomerization reactions (Transfer of a primary H-atom)
	1-4 H-transfer	1-5 H-transfer	1-6 H-transfer
Primary radical	10^11exp(-20,600/RT)	1.58×10^10exp(-14,500/RT)	3.16×10^9exp(-19,500/RT)
Secondary radical	10^11exp(-21,600/RT)	1.58×10^10exp(-15,500/RT)	3.16×10^9exp(-20,500/RT)
Tertiary radical	10^11exp(-22,100/RT)	1.58×10^10exp(-16,000/RT)	3.16×10^9exp(-21,000/RT)
Alkyl radical decomposition reactions (to form primary radicals)
	Primary radical	Secondary radical	Tertiary radical
	10^14exp(-30,000/RT)	10^14exp(-31,000/RT)	10^14exp(-31,500/RT)
Corrections of decomposition rates to form
Methyl radical	Secondary radical	Tertiary radical	Alkyl radical
exp(-2,500/RT)	exp(1,500/RT)	exp(2,500/RT)	0.316×exp(8,000/RT)

Table 1.

Reference kinetic parameters of pyrolysis reactions [26].

Figure 4.

Flowchart of automatic generation and reduction of reaction network.

$$\frac{d\left[c\right]}{dt}=r_{formation}-r_{consumption}\approx 0$$E7

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

$$r_j+{\displaystyle \sum _{i\ne j}{k}_{i,j}{R}_i}-\left|{\displaystyle \sum _{i\ne j}{k}_{j,i}+k_{d_i}}\right|R_j=0$$E8

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 $\to$ 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.

1. Molecular-level model

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.

Reaction	A	n	E; (kcal/mol)	Reaction	A	n	E; (kcal/mol)
CH3J+CH3J $\leftrightarrow$ C2H6	8.26E+17	-1.4	1	BUT-2+C2H5J $\leftrightarrow$ C4H7J+C2H6	3.36E+12	0	12.39
CH3J+C2H5J $\leftrightarrow$ C3H8	1.37E+13	0	0	BUT-1+HJ $\leftrightarrow$ C4H7J+H2	1.40E+04	2.36	1.11
C2H5J+C2H5J $\leftrightarrow$ n-C4H10	1.15E+13	0	0	BUT-2+HJ $\leftrightarrow$ C4H7J+H2	2.60E+06	2.38	2.8
C2H5J+HJ $\leftrightarrow$ C2H6	1.00E+14	0	0	BUT-1+C3H7J $\leftrightarrow$ C4H7J+C3H8	1.94E+02	2.96	6.79
C2H4+HJ $\leftrightarrow$ C2H5J	1.55E+14	0	2.8	BUT-2+C3H7J $\leftrightarrow$ C4H7J+C3H8	1.72E+12	0	12.29
C2H6+HJ $\leftrightarrow$ C2H5J+H2	2.27E+08	1.75	7.51	BUT-1+C3H7J $\leftrightarrow$ C4H7J+C3H8	3.12E-04	4.31	3.39
HJ+HJ $\leftrightarrow$ H2	1.09E+11	0	1.5	BUT-2+C3H7J $\leftrightarrow$ C4H7J+C3H8	3.36E+12	0	12.39
C3H7J+HJ $\leftrightarrow$ C3H8	2.00E+13	0	0	BUT-1+C2H3J $\leftrightarrow$ C4H7J+C2H4	5.11E+00	3.59	5.06
CH3J+C3H7J $\leftrightarrow$ i-C4H10	6.64E+14	-0.57	0	BUT-2+C2H3J $\leftrightarrow$ C4H7J+C2H4	4.64E+13	0	7.5
C3H8+C2H5J $\leftrightarrow$ C3H7J+C2H6	3.08E+03	2.66	10.1	BUT-1+C4H7J $\leftrightarrow$ C4H7J+BUT-1	3.12E-04	4.31	3.39
C3H8+HJ $\leftrightarrow$ C3H7J+H2	1.60E+08	1.69	4.78	BUT-2+C4H7J $\leftrightarrow$ C4H7J+BUT-1	3.36E+12	0	12.39
C3H6+HJ $\leftrightarrow$ C3H7J	2.01E+13	0	2.1	C4H7J $\leftrightarrow$ C4H7J	2.82E+08	1.28	27.89
C3H7J+HJ $\leftrightarrow$ C3H8	1.00E+14	0	0	C4H7J+BUT-1 $\leftrightarrow$ BUT-2+C4H7J	3.12E-04	4.31	3.39
CH3J+C3H7J $\leftrightarrow$ n-C4H10	3.37E+13	0	0	C2H3J+C2H3J $\leftrightarrow$ C4H6	7.23E+13	0	0
C2H4+CH3J $\leftrightarrow$ C3H7J	1.94E+13	0	9	C4H6+HJ $\leftrightarrow$ C4H7J	5.70E+13	0	4.3
C3H8+C2H5J $\leftrightarrow$ C3H7J+C2H6	3.95E+03	2.71	12.91	C4H6+HJ $\leftrightarrow$ C4H7J	7.42E+13	0	1
C3H8+HJ $\leftrightarrow$ C3H7J+H2	3.77E+08	1.75	7.51	C3H5J+HJ $\leftrightarrow$ C3H6	5.00E+13	0	0
C3H6+HJ $\leftrightarrow$ C3H7J	1.18E+13	0	3.8	C3H5J+CH3J $\leftrightarrow$ BUT-1	1.00E+13	0	0
C3H7J+C3H8 $\leftrightarrow$ C3H8+C3H7J	3.08E+03	2.66	10.1	C3H6+CH3J $\leftrightarrow$ C3H5J+CH4	2.41E+02	2.92	7.16
C3H7J $\leftrightarrow$ C3H7J	1.94E+10	0.89	35.78	C3H6+C2H5J $\leftrightarrow$ C3H5J+C2H6	1.01E-04	4.75	4.13
C2H6+CH3J $\leftrightarrow$ C2H5J+CH4	1.67E+06	1.9	11.04	C3H6+HJ $\leftrightarrow$ C3H5J+H2	1.30E+06	2.38	2.8
C3H8+CH3J $\leftrightarrow$ C3H7J+CH4	1.67E+06	1.9	11.04	C3H6+C3H7J $\leftrightarrow$ C3H5J+C3H8	8.61E+11	0	12.29
C3H8+CH3 $\leftrightarrow$ C3H7J+CH4	2.90E+06	1.77	8.53	C3H6+C3H7J $\leftrightarrow$ C3H5J+C3H8	1.01E-04	4.75	4.13
HJ+CH3J $\leftrightarrow$ CH4	1.93E+14	0	0.27	C3H6+C2H3J $\leftrightarrow$ C3H5J+C2H4	2.32E+13	0	7.5
HJ+CH4 $\leftrightarrow$ H2+CH3J	1.22E+08	1.87	10.59	C3H6+C4H7 $\leftrightarrow$ C3H5J+BUT-1	1.01E-04	4.75	4.13
C2H3J+HJ $\leftrightarrow$ C2H4	5.36E+14	0	0.98	C3H5J+BUT-1 $\leftrightarrow$ C3H6+C4H7J	3.12E-04	4.31	3.39
C2H3J+CH3J $\leftrightarrow$ C3H6	7.23E+13	0	0	C3H5J+BUT-2 $\leftrightarrow$ C3H6+C4H7J	3.47E+13	0	21.39
C2H3J+C2H5J $\leftrightarrow$ BUT-1	1.00E+13	0	0	C4H7J+C2H3J $\leftrightarrow$ C6H10	1.00E+13	0	0
C2H3J+CH4 $\leftrightarrow$ C2H4+CH3J	5.36E+14	0	13.69	C3H5J+C3H5J $\leftrightarrow$ C6H10	1.00E+13	0	0
C2H3J+C2H6 $\leftrightarrow$ C2H4+C2H5J	5.30E+04	2.44	5.5	C4H6+CH3J $\leftrightarrow$ C5H9J	1.46E+13	0	5.7
C2H3J+H2 $\leftrightarrow$ C2H4+HJ	8.98E+12	0	10.3	C2H4+C3H5J $\leftrightarrow$ C5H9J	1.43E+12	0	14.39
C2H3J+C3H8 $\leftrightarrow$ C2H4+C3H7J	1.02E+03	3.1	8.82	C5H9 $\leftrightarrow$ C5H9J	3.32E+08	1.1	29.39
C2H3J+C3H8 $\leftrightarrow$ C2H4+C3H7J	9.30E+04	2.44	5.5	C4H9J+H $\leftrightarrow$ n-C4H10	1.00E+14	0	0
C4H7J+HJ $\leftrightarrow$ BUT-1	1.00E+14	0	0	n-C4H10+CH3 $\leftrightarrow$ C4H9J+CH4	1.67E+06	1.9	11.04
BUT-1+CH3J $\leftrightarrow$ C4H7J+CH4	8.34E+05	1.9	11.04	C4H9J+C2H6 $\leftrightarrow$ n-C4H10+C2H5J	1.93E-05	5.28	7.78
C4H7J+C2H6 $\leftrightarrow$ BUT-1+C2H5J	1.93E-05	5.28	7.78	C2H4+C2H5J $\leftrightarrow$ C4H9J	3.98E+03	2.44	5.37
BUT-1+HJ $\leftrightarrow$ C4H7J+H2	1.88E+08	1.75	7.51	n-C4H10+HJ $\leftrightarrow$ C4H9J+H2	3.77E+08	1.75	7.51
C4H7J+C3H8 $\leftrightarrow$ BUT-1+C3H7J	3.08E+03	2.66	10.1	BUT-1+HJ $\leftrightarrow$ C4H9J	1.18E+13	0	3.8
C4H7J+C3H8 $\leftrightarrow$ BUT-1+C3H7J	3.95E+03	2.71	12.91	C4H9J+C3H8 $\leftrightarrow$ n-C4H10+C3H7J	3.08E+03	2.66	10.1
BUT-1+C2H3J $\leftrightarrow$ C4H7J+C2H4	4.65E+04	2.44	5.5	C4H9J+C3H8 $\leftrightarrow$ n-C4H10+C3H7J	3.95E+03	2.71	12.91
C2H4+C2H3J $\leftrightarrow$ C4H7J	4.18E+10	0	5.2	n-C4H10+C2H3 $\leftrightarrow$ C4H9J+C2H4	9.30E+04	2.44	5.5
C4H7J+HJ $\leftrightarrow$ BUT-1	5.00E+13	0	0	n-C4H10+C4H7J $\leftrightarrow$ C4H9J+BUT-1	3.95E+03	2.71	12.9
C4H7J+HJ $\leftrightarrow$ BUT-2	5.00E+13	0	0	C4H9J+BUT-1 $\leftrightarrow$ n-C4H10+C4H7J	3.12E-04	4.31	3.39
BUT-1+CH3J $\leftrightarrow$ C4H7J+CH4	5.46E+13	0	10.39	C4H9J+BUT-2 $\leftrightarrow$ n-C4H10+C4H7J	3.36E+12	0	12.39
BUT-2+CH3J $\leftrightarrow$ C4H7J+CH4	4.82E+02	2.92	7.16	C4H9J+C3H6 $\leftrightarrow$ n-C4H10+C3H5J	1.01E-04	4.75	4.13
BUT-1+C2H5J $\leftrightarrow$ C4H7J+C2H6	3.12E-04	4.31	3.39

Table 2.

Elementary reaction network for light hydrocarbon steam cracking.

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).

$$r_c={\displaystyle \sum }_{i=1}^9{A}_i\exp \left(-\frac{E_i}{RT}\right)c_i$$E9

In eq. (9): 1=C₂H₂; 2=C₂H₄; 3=C₃H₆; 4=1-C₄H₈; 5=C₄H₆; 6=C₆H₆; 7=C₇H₈; 8=xylene; 9=C₂H₃-C₆H₅.

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.

No.		1	2	3	4	5	6	7	8	9	10	11	12	13
Feed
D15		0.7091	0.7106	0.7109	0.7099	0.7122	0.694	0.7099	0.6981	0.7109	0.6796	0.7121	0.7115	0.7109
M		94.1	94.76	94.68	94.29	94.31	94.43	94.6	93.67	94.74	95.26	95.26	95.41	94.83
T	°C	99.3	101.7	101.6	103.38	101.7	99.6	96.4	99	100.9	101.2	103.3	98.2	97
K		12.45	12.45	12.45	12.48	12.41	12.71	11.17	12.64	12.44	13.02	12.44	12.4	12.39
Operating conditions
F	t/h	6.3686	5.413	5.4166	6.0792	5.1948	5.2302	4.9632	4.9812	5.3286	5.397	5.3726	5.431	5.5722
W/O	t/t	0.5	0.53	0.53	0.5	0.53	0.53	0.54	0.53	0.52	0.54	0.53	0.52	0.53
CIT	°C	639.6	639.4	639.5	633.4	640.9	639.8	642.5	641.4	636.5	634.8	638.8	640.8	639.1
COT	°C	871.6	867.6	867.2	872.6	865.6	869.2	865.8	861.6	865.4	863.4	868.4	865.8	867
CIP	Mpa	0.366	0.314	0.312	0.345	0.641	0.303	0.292	0.289	0.302	0.312	0.311	0.313	0.326
COP	Mpa	0.044	0.043	0.042	0.039	0.043	0.048	0.044	0.04	0.044	0.044	0.043	0.042	0.045
Yields of major products
H2	wt%	1.07	0.99	1.09	1.10	0.99	1.07	0.87	1.01	0.93	1.15	1.19	1.06	1.15
CH4	wt%	16.55	15.62	16.73	16.83	15.39	16.73	13.90	15.76	15.15	18.09	18.80	16.98	17.64
C2H4	wt%	30.65	28.62	29.59	30.03	28.57	30.52	29.23	28.45	30.91	30.64	30.82	30.55	29.56
C2H6	wt%	3.79	3.44	3.54	3.56	3.34	3.64	3.49	3.69	3.61	3.63	3.65	3.68	3.60
C3H6	wt%	16.13	17.38	16.51	16.09	14.88	14.95	15.72	16.95	17.22	15.89	16.50	17.08	16.72
C3H8	wt%	0.90	1.74	0.87	0.82	0.83	0.69	0.99	0.88	1.53	0.96	1.05	1.21	1.17
C4H6	wt%	6.02	6.13	6.09	6.03	5.74	5.76	5.51	6.06	5.63	5.62	5.80	5.57	5.48
N-C4H8	wt%	1.41	1.59	1.70	1.54	1.40	1.31	1.46	2.07	1.38	1.37	1.43	1.58	1.57
I-C4H8	wt%	2.66	2.83	2.98	2.70	2.38	2.32	2.57	3.18	2.79	2.67	2.88	3.09	3.09
others	wt%	20.82	21.66	20.9	21.3	26.48	23.01	26.26	21.95	20.85	19.98	17.88	19.2	20.02

Table 3.

Industrial data for feedstock, operating conditions and yields of major products.

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₁, COT₂,..., COT_n].

The optimization model is summarized below.

$$\max \frac{1}{N}{\displaystyle \sum }_{j=1}^M{\displaystyle \sum }_{i=1}^N{\omega }_t{y}_{j,t}^o$$E10

Subject to:

$$\begin{array}{ccc}\frac{d{y}_{m,t}}{dL}=\frac{S_t}{V_t}{\displaystyle \sum _i{\nu }_{i,m}{r}_{i,t}}& \forall m\in \left\{1,2,\cdots ,M\right\}& \forall t\in \left\{1,2,\cdots ,N\right\}\end{array}$$E11

$$\begin{array}{cc}\frac{d{P}_t}{dL}=-\frac{f\cdot E\left(L\right)\cdot G_t^2}{5.07\times {\rho }_t\cdot D_{i,t}\times {10}^4}& \forall t\in \left\{1,2,\cdots ,N\right\}\end{array}$$E12

$$\begin{array}{cc}\frac{d{T}_t}{dL}=\frac{k_t\pi D_o\left(T_{w,t}-T_t\right)-{\displaystyle \sum _m\Delta H_{fm,t}^0\cdot \frac{d{y}_{m,t}}{dL}}}{{\displaystyle \sum _m{c}_{pm,t}{y}_{m,t}}}& \forall t\in \left\{1,2,\cdots ,N\right\}\end{array}$$E13

$$\begin{array}{cc}{D}_{i,t+1}=D_{i,t}-2\Delta {\delta }_{c,t}=D_{i,t}-2\cdot \frac{r_{c,t}}{{\rho }_{c,t}}& \forall t\in \left\{1,2,\cdots ,N-1\right\}\end{array}$$E14

$$\begin{array}{cc}{r}_{c,t}={\displaystyle \sum }_{i=1}^9{A}_i\exp \left(-\frac{E_i}{RT}\right)c_{i,t}& \forall t\in \left\{1,2,\cdots ,N\right\}\end{array}$$E15

$$\begin{array}{cc}\underset{t}{\max }{T}_{w,t}\le T_M& \forall t\in \left\{1,2,\cdots ,N\right\}\end{array}$$E16

$$\begin{array}{cc}\left|CO{T}_{t+1}-CO{T}_t\right|\le \delta & \forall t\in \left\{1,2,\cdots ,N-1\right\}\end{array}$$E17

$$\begin{array}{cc}LB\le CO{T}_t\le UB& \forall t\in \left\{1,2,\cdots ,N\right\}\end{array}$$E18

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.

Figure 5.

Simulation results compared with industrial data.

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.

$$d_t=f\left(\text{feed},\text{period }\text{}\text{number},CO{T}_1,CO{T}_2,\cdots ,CO{T}_t\right)$$E19

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).

$$\begin{array}{l}\begin{array}{cc}\delta d_t=d_{t+1}-d_t=f\left(t,CO{T}_t\right)=a\times CO{T}_t^2+b\times t^2+c\times CO{T}_t+d\times CO{T}_t+e\times t+f& \end{array}\\ \forall t\in \left\{1,2,\cdots ,N-1\right\}\end{array}$$E20

$$d_t={\displaystyle \sum _{t=1}^i\delta d_t}={\displaystyle \sum _{t=1}^i{f}_t\left(t,CO{T}_t\right)}$$E21

Figure 6.

Optimization procedure with multi-period steam cracking simulation.

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.

Figure 7.

Coke thickness from the surrogate model and original multi-period model.

Figure 8.

Optimization using a parallel simulation procedure.

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.

Figure 9.

Detailed scheme for parallel simulation.

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.

Figure 10.

Comparison between a warm start for serial and parallel simulation, based on the multi-period steam cracking model.

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.

Feed
Molecular weight	95.26
H/C	2.11
	P		I			N			A
PINA	0.3096		0.3223			0.2856			0.0825
	IBP	10%		30%	50%		70%	90%		EBP
ASTM (°C)	41.2	58.2		80.1	99.0		119.6	148.9		175.9
Operating condition
Feed rate (t/h)	5.4
Water/Oil (t/t)	0.54
Coil inlet temperature (°C)	636.5
Coil inlet pressure (Mpa)	0.31
Coil outlet pressure (Mpa)	0.04
Optimization variables
Objective function	$\max \frac{1}{10}\left(0.52{\displaystyle \sum _{t=1}^{10}{y}_{C_2{H}_4,t}^o}+0.48{\displaystyle \sum _{t=1}^{10}{y}_{C_3{H}_6,t}^o}\right)$
Furnace runtime (day)	50
NP	8
TM (°C)	1000.0
δ (°C)	5.0
LB (°C)	855.0
UB (°C)	872.0

Table 4.

Supporting information for the steam cracking optimization model.

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),

$$\begin{array}{cc}\delta d_t=d_{t+1}-d_t=f\left(t,CO{T}_t\right)& \forall t\in \left\{1,2,\cdots ,N-1\right\}\end{array}$$E22

Figure 11.

COT optimization results for serial and parallel simulation.

Figure 12.

C₂H₄ and C₃H₆ optimization results for serial and parallel calculation.

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

$$\begin{array}{cc}\delta d_t=d_{t+1}-d_t=f\left(t,CO{T}_1,CO{T}_2,\cdots ,CO{T}_t\right)& \forall t\in \left\{1,2,\cdots ,N-1\right\}\end{array}$$E23

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.

	COT			C2H4			C3H6
Period	Serial (°C)	Parallel (°C)	Error (%)	Serial (wt%)	Parallel (wt%)	Error (%)	Serial (wt%)	Parallel (wt%)	Error (%)
1	872.00	872.00	0.00	32.13	32.13	0.01	16.54	16.54	0.00
2	872.00	872.00	0.00	32.22	32.24	0.05	16.55	16.55	0.00
3	872.00	872.00	0.00	32.31	32.36	0.14	16.57	16.57	0.01
4	872.00	872.00	0.00	32.42	32.48	0.21	16.59	16.59	0.02
5	872.00	872.00	0.00	32.56	32.65	0.27	16.58	16.56	0.14
6	872.00	872.00	0.00	32.70	32.84	0.44	16.55	16.49	0.34
7	870.58	871.38	0.09	32.84	33.02	0.54	16.46	16.40	0.34
8	865.58	866.38	0.09	32.89	33.07	0.54	16.31	16.22	0.52
9	860.58	861.38	0.09	32.98	33.09	0.35	16.11	16.04	0.39
10	855.58	856.38	0.09	33.11	33.09	0.04	15.87	15.88	0.07
Average error	-	-	0.04	-	-	0.26	-	-	0.18

Table 5.

Summarized optimization results for serial and parallel calculation.

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.

Acknowledgments

The authors gratefully acknowledge the financial support of the National Basic Research Programme (No. 2012CB720500) and the NSFC (No. 21306100).