Ind. Eng. Chem. Res. 2009, 48, 9523–9533
9523
Multiobjective Optimization of an Industrial LPG Thermal Cracker using a First Principles Model S. R. Nabavi,†,‡ G. P. Rangaiah,*,† A. Niaei,‡ and D. Salari‡ Department of Chemical and Biomolecular Engineering, National UniVersity of Singapore, Engineering DriVe 4, Singapore 117576, Republic of Singapore, and Department of Applied Chemistry and Chemical Engineering, UniVersity of Tabriz, 51666-16471, Tabriz, Iran
Thermal cracking of hydrocarbons in the presence of steam is the most important process for the production of ethylene and propylene, which are the raw materials for many petrochemicals. In this work, multiobjective optimization (MOO) of an industrial LPG (liquefied petroleum gas containing mainly propane and butane) thermal cracker is studied using the elitist nondominated sorting genetic algorithm adapted with the jumping gene operator, NSGA-II-aJG. A first principles model based on free radical mechanism for LPG thermal cracking is employed for optimization. Several bi- and triobjective optimization problems are solved; these problems involved maximization of annual ethylene/propylene production, selectivity and run length, and minimization of severity and total heat duty per year. Feed flow rate, steam ratio, inlet temperature, coil outlet temperature, and pressure are the decision variables. MOO of the LPG thermal cracker provides a range of optimal operating conditions and objective values, with the triobjective optimization giving more optimal solutions compared to those from the biobjective optimization for the objectives considered. With this detailed knowledge, a suitable operating point can be selected based on specific requirements in the particular plant. 1. Introduction Steam cracking is the most important technology for the production of light olefins and aromatics. Hydrocarbons such as ethane, propane, liquefied petroleum gas (LPG), virgin naphtha, and gas oil mixed with steam are fed into tubular reactors (cracking coils) at high temperature.1 They are cracked into a combination of olefins, aromatics, pyrolysis fuel oil, and other heavier hydrocarbons. The reactions take place in several cracking coils placed vertically inside the fire-box and heated by gas-fueled (hydrogen and methane mixture) burners. The products from the reaction are quenched in a transfer line exchanger (TLE). Coke is formed on both the cracking coil inner surface and in the TLE. It reduces the olefins selectivity mainly because of the increased pressure drop and hence pressure in the cracking coils. The coke formation also leads to increased tube skin temperature and to the operational constraint (namely, maximum limit on tube skin temperature) for shutdown of the furnace. Coke is burned off with steam and air during 1-4 days. This decoking operation is required every 40-80 days depending on furnace design and operation approach. Ethylene, the main product from steam cracking is a very reactive intermediate, which can undergo all typical reactions of a short-chain olefin. Hence, ethylene has gained importance as a chemical building block. Most of the ethylene is used to produce ethylene oxide, ethylene dichloride, and polyethylene. Significant amounts of ethylene are also used to make ethyl benzene, oligomer products, acetaldehyde, and vinyl acetate.2 The most valuable byproduct in ethylene production is propylene, the second most important feedstock for petrochemical products. The major propylene derivatives include polypropylene, acrylonitrile, propylene oxide, cumene/phenol, oxo alcohols, acrylic acid, isopropyl alcohol, oligomers, and other * To whom correspondence should be addressed. E-mail: chegpr@ nus.edu.sg. † National University of Singapore. ‡ University of Tabriz.
intermediates (www.plastemart.com, Technical Articles & Reports on Plastic Industries, accessed September 2008). Therefore, any improvement in the production of ethylene and propylene can generate significant revenue for petrochemical industries. There have been a number of studies on the optimization of thermal cracking plants. Optimal operation of an ethylene plant at variable feed conditions has been reported by Eliceche et al.3 Optimization of the total plant including furnace and downstream units was studied for variable feed flow rate and composition. Nonlinear programming (NLP) was used for maximization of the profit. Active constraints related to two pieces of equipmentsthe furnace and ethylene-ethane splitters indicated plant bottlenecks. Dynamic optimization of the production period of thermal cracking with consideration of coke formation in cracking coil and TLE was carried out by Edwin and Balchen.4 Optimal time-dependent trajectories of feed rate, steam to hydrocarbon ratio, and conversion were found for maximization of net profit. Discretization of the trajectories and sequential quadratic programming (SQP) were employed to solve the optimization problem.4 Dynamic programming was applied to the optimal control of propane thermal cracker.5,6 The molecular reaction scheme was used for modeling propane cracking in the industrial reactor, and the optimal temperature profile for profit maximization was obtained. Schulz et al.7 formulated scheduling of parallel ethanefed cracking furnace shutdowns in an ethylene plant as a multiperiod mixed-integer NLP (MINLP) with discrete time representation. Nonlinear correlations mainly for furnace operation and separation train were included for the entire process at each time interval. The objective was profit function based on revenue (from sale of ethylene, propylene, propane, butane, gasoline, and residue gas) and cost (of feed, natural gas for boilers and furnaces, inventory, and cleanup of furnaces). In all the above cases, optimization was carried out for single objective even though multiple objectives are possible. For example, simultaneous maximization of some objectives such as feed conversion and olefin yield is useful. Solutions of MOO
10.1021/ie801409m CCC: $40.75 2009 American Chemical Society Published on Web 03/03/2009
9524
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
are conceptually different form single-objective optimization. In MOO, there may not be a single solution that is the best with respect to all objectives. Instead, there may be a set of optimal solutions that are equally good, which are known as Pareto-optimal solutions. This set of solutions is such that when we go from any one solution to another in the set, at least one objective improves and at least one other worsens. Bhaskar et al.8 reviewed MOO background, different methods, and their applications in chemical engineering. Applications of nondominated sorting genetic algorithm (NSGA) and its variants in chemical reaction engineering were reviewed by Nandasana et al.9 The latest review of Masuduzzaman and Rangaiah10 shows nearly 100 chemical engineering applications of MOO reported since the year 2000. However, there have been only a few studies on MOO of thermal crackers. One of them is by Tarafder et al.11 for an industrial ethylene reactor using ethane feed. Simulation of the reactor was based on a free radical mechanism without coke formation reactions. The elitist NSGA (NSGA-II) was used for optimization of the cracking reactor. Different MOO cases involving maximization of ethane conversion, ethylene selectivity, and flow rate, as well as design and operation stages, were considered, and the effect of both design and operation variables on the optimal solutions was discussed.11 Recently, Li et al.12 investigated MOO of a naphtha thermal cracker for simultaneous maximization of ethylene and propylene yield, using the multiobjective particle swarm optimization. A neural network model trained on industrial data was employed for modeling the reactor; this model predicts the yield of ethylene and propylene for given coil outlet temperature, coil outlet pressure, and steam-to-naphtha ratio. More recently, Gao et al.13 optimized the periodic operation of a naphtha cracker for maximizing ethylene and propylene yield using a parallel hybrid algorithm based on NSGA-II and SQP. The kinetic model was based on a molecular mechanism, and an operating period of 20 days was assumed for optimization. Motivated by the limited number of studies on MOO of the industrially important thermal crackers, the main aim of this study is to optimize a thermal cracker with liquefied petroleum gas (LPG) as the feed, for multiple objectives using a detailed, first principles model. Unlike in the previous works, the first principles model is based on a free radical mechanism with coke formation reactions for realistic simulation of thermal crackers. Further, simulation is carried out for one complete run from the start to decoking time. The elitist nondominated sorting genetic algorithm adapted with the jumping gene operator, NSGA-II-aJG,14 is employed for MOO. Several bi- and triobjective cases involving important and conflicting objectives were considered for MOO of a LPG thermal cracker. The results are presented and discussed to highlight their validity and the salient features. The rest of this paper is organized as follows. In section 2, reaction mechanisms and modeling of thermal crackers for gaseous feeds are briefly reviewed. The simulation and optimization procedures adopted for this study on an LPG thermal cracker are described in section 3. Formulation of MOO problems in section 4 covers selection of objectives, decision variables, and constraints in the optimization problems studied. In section 5, results from the optimization of several combinations of two and three objectives are presented and discussed. Finally, conclusions of this study are given in section 6.
2. Reaction Mechanism and Modeling of the Cracking Process For realistic simulation of the thermal cracking reactors, a complex kinetic model featuring coke formation needs to be combined with models for heat transfer, fire-box profiles, and fluid dynamics in the furnace. The reaction mechanism, the central part of the simulation, can be overall, molecular, and free radical mechanisms, of which the last is the most detailed and perhaps the most accurate. Several free radical and molecular mechanisms have been reported for light hydrocarbon cracking in the open literature, and some of these are briefly reviewed here. A molecular mechanism approximating the free radical nature of ethane thermal cracking was proposed by Froment et al.15 and Sundaram and Froment.16 Molecular mechanisms lead to nonstiff differential equations that are easier to solve whereas free radical mechanisms lead to stiff differential equations that are difficult to solve. Sundaram and Froment17 developed a free radical mechanism for ethane cracking involving 49 reactions. In this scheme, the reaction mechanism was simplified by lumping products heavier than C5H10 with very small yields, into a single component C5+. Thermal cracking of pure propane and its mixtures with other hydrocarbons has also been studied. A free radical mechanism with eight reactions but without coke formation was developed by Herriote et al.,18 based on data from an experimental tubular reactor. Van Damme et al.19 developed a molecular scheme with 16 elementary reactions for thermal cracking of propane and propane-propylene mixture. The reaction mechanism was validated against the data from a pilot plant and then used for simulating an industrial reactor. Sundaram and Froment20 proposed a molecular model including coke formation for propane cracking, based on data from a laboratory reactor, and then used it for simulating an industrial reactor. In most of the above cases, the kinetic model was inserted into the reactor model and the simulation was carried out for predicting product distribution, temperature and pressure profiles, and run length of industrial cracking furnace (i.e., without considering fire-box modeling and simulation). Plehiers et al.21 successfully simulated the ethane cracking furnace using a combined fire-box and reactor model including coke formation reactions. The simulation approach provided temperature and heat flux distribution for each increment in run time. Pseudosteady state was assumed for coke formation; i.e., the rate of coke formation is constant in each time interval, and a time interval of 65 h was found to be suitable. Coke formation is an unwanted process that takes place during thermal cracking. It reduces heat transfer from the furnace to the process gas and also reduces the cross section of the reactor. All these dictate the furnace run length and shutdown for reactor decoking, which involves a significant loss in production capacity and should be avoided. LPG is a mixture of mainly propane, isobutane, and n-butane. A kinetic model of LPG thermal cracking including coke formation reactions was developed by Towfighi et al.22 The free radical mechanism developed contains 146 elementary reactions with molecular and free radical species. This reaction scheme was inserted into the reactor model, and the simulation of a cracking furnace for an entire run length was successfully carried out. An independently simulated, external heat flux profile was used in the reactor model. By using this approach, the computational time was reduced and manipulation of the heat flux profile can be carried out easily. The present study on MOO of the LPG thermal cracker is based on this first principles model
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
9525
Table 1. Values of NSGA-II-aJG Parameters Used in This Study
Figure 1. Flowchart for reactor simulation.
of LPG thermal cracking and the approach of independent simulation of the heat flux profile of Towfighi et al.22 3. Simulation and Optimization Procedures For the simulation of an LPG cracking reactor, the governing mass, energy, and momentum balance equations should be solved simultaneously. These equations and the continuity equation for coke deposition are listed in Appendix A. Owing to the significant stiffness of these equations arising from the large differences in concentration gradients between free radicals and molecular species, Gear’s method was used for solving the differential equations. A flowchart of the iterative scheme for thermal cracker simulation is given in Figure 1. Predicting the run length of the reactor requires the simultaneous solution of numerous differential equations (eqs A-1, A-4, A-5, and A-8 in Appendix A). By assuming pseudosteady state for cracking reactions since they are much faster than coke formation reactions, the run time of the cracking furnace can be incremented in a stepwise manner. During the simulations, the coil outlet temperature (COT) and pressure (COP) should be kept constant at the specified values as per the industrial practice (or by the optimizer during the optimization). Therefore, two iterative loops are needed for maintaining these values (see Figure 1). Owing to coke formation on tube/coil internal surface and the poor thermal
parameter
value
number of chromosomes in the population, Npop number of generations, Ngen length of a substring, lsustr total length (number of binaries) of a chromosome, lchrom length of the replacing jumping gene, laJG crossover probability, Pc mutation probability, Pm jumping probability for the JG operator, PaJG seed for random number generator, Sr
50 100 30 150 40 0.8 0.01 0.8 0.6
conductivity of the coke layer, the process gas temperature decreases. To maintain the COT at its set point, the heat input into the reactor should be increased with time. Hence, the heat flux profile is changed to achieve the specified COT. Furthermore, coke formation decreases the internal cross section of tubes and results in a higher pressure drop over the reactor. To maintain the required COP, the inlet pressure has to be increased with time. Finally, the furnace has to be stopped for decoking when the external skin temperature of tubes reaches the maximum value (namely, 1100 °C) imposed by the tube material (Figure 1). A time increment of 48 h of furnace run time was found to be sufficiently small to justify pseudosteady state.22 The kinetic model of Towfighi et al.22 was validated using the pilot plant data, which was designed and assembled by the Olefin Research Group. Details of this pilot plant are available at http:// www.Modares.ac.ir/English/faculties/eng/olefin/index.htm. Towfighi et al.22 used their model along with independently simulated heat flux profile, for simulating an industrial LPG cracker in the olefin plant of Abadan Petrochemicals, Iran, for different feed composition and severity index; the predicted yields of products were comparable to the plant data. Many methods are available for solving MOO problems. Of these, methods based on genetic algorithms have found many applications in chemical engineering.10 Several extensions of simple genetic algorithm have been developed to solve MOO problems.23 One of the popular algorithms is the elitist NSGAII.24 It was modified to NSGA-II-JG by Kasat and Gupta25 by including a modified mutation operator based on the concept of jumping genes (JG) in natural genetics. This is a macromacro mutation and counteracts the decrease in the diversity created by elitism. Kasat and Gupta25 found that the CPU time required for MOO is reduced 5-fold when NSGA-II-JG was used as compared to when NSGA-II was used. Guria et al.14 used another JG adaptation of NSGA-II (referred to as NSGA-II-aJG) for MOO of reverse osmosis desalination unit. They applied different versions of NSGA-II and concluded that NSGA-II-aJG converges faster than others tried. Agrawal et al.26 applied binary-coded NSGA-II-JG and NSGA-II-aJG for MOO of an industrial low density polyethylene reactor at operation stage. Later, Agrawal et al.27 studied MOO of the same reactor at design stage. In both these studies, it was found that the binary-coded NSGA-II-aJG and NSGA-II-JG performed better than the binary-coded NSGA-II near the hard end-point constraints. Hence, NSGA-II-aJG was chosen for MOO in this work. The working procedure and flowchart of NSGA-II-aJG can be found in ref 26. Similar to other versions of NSGA, NSGA-II-aJG has some parameters that should be selected for optimization. In this study, values of these parameters were based on those used by Agrawal et al.26,27 and several trials with different values of parameters. Finally, the parameter values in Table 1 were used for all optimization results presented in the next section.
9526
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
4. Multiobjective Optimization In most optimization studies on thermal cracking and other processes, the common objective is profit or costs. However, the profit and costs encompass several factors such as demand and prices that fluctuate with time and location, downstream processing costs, etc. Hence, this study considers different kinds of objectives that increase the scope of making profit rather than maximizing the profit itself. They are also fundamental in characterizing the cracker performance. The objectives considered for LPG cracker optimization are maximization of ethylene and propylene production per year, selectivity index, and run length and minimization of severity index and total heat duty per year. Annual ethylene and propylene production and heat duty are obviously important in many industries. Run length in steam cracking is important to avoid frequent disruptions to the plant. The ethylene selectivity index is the ratio of ethylene weight fraction to ethane weight fraction, and the severity index is the ratio of some lighter fraction (such as propane, propylene, propadiene, ethane, ethylene, acetylene, methane, and hydrogen) to the propylene fraction. Van Geem et al.28 concluded that, for similar feed stocks, the product yields from different reactor coil geometry can be characterized by these two indices. Hence, the objectivessselectivity and severity indicessare not dependent on reactor configuration, and the optimization results are applicable to cracking of similar LPG feed stocks in different coil configurations. We considered several combinations of two and three objectives from the above seven objective functions. In this paper, we present and discuss results for five cases, each with two or three objectives. These were selected for the conflicting nature of objectives and for coherent discussion of results from bi- and triobjective optimization. The five cases to be presented in sections 5.1 and 5.2 are the following: • Case I: maximization of ethylene and propylene production per year • Case II: maximization of ethylene selectivity index and minimization of severity index • Case III: maximization of ethylene production and minimization of heat duty per year • Case A: maximization of ethylene selectivity index and run length and minimization of severity index • Case B: maximization of ethylene and propylene production and minimization of heat duty As the selected optimization routine, NSGA-II-aJG, is based on solving a maximization problem, the minimization objective function must be modified using a suitable transformation. For example, a minimization objective, J, can be transformed to a maximization objective, I, by I)
1 J+1
(1)
The variables, which significantly affect the thermal cracking reactor performance and can be adjusted in an industrial reactor system, were chosen as the decision variables. These variables and their ranges are as follows. 10 e Fin e 14 t/(h furnace)
(2)
800 e COT e 850 °C
(3)
1.5 e COP e 2.5 bar
(4)
0.2 e SR e 0.4 kg steam/kg feed
(5)
500 e Tin e 700 °C
(6)
Table 2. Basic Data of the Cracking Coil and Fire-Box of an LPG Cracker Furnace fire-box height (m) length (m) depth (m) no. of burners
cracking coil (reactor) 13.3 11.7 3.0 112
fuel composition (% mol) H2 CH4 C2H6 excess air ) 15%
3.6 95.4 1
length (m) internal diameter (mm) external diameter (mm)
98 128 141
feed composition (wt %) C2H6 C3H8 nC4 iC4
5.18 16.9 59.6 18.32
The feed flow rate is for each furnace with four coils, and so, the inlet flow rate (Fin) is divided equally to these four coils. Basic data of the fire-box, cracking coil, and feed and fuel compositions are summarized in Table 2. COT is a key variable in the thermal cracking operation; a high value of COT results in overcracking, which can lead to premature furnace shutdown due to excessive coke formation inside the reactor. Yet, undercracking of the feed due to low COT reduces overall production and hence there exists an optimal COT. The lower and upper bounds on COT in eq 3 are based on typical industrial practice. The COP is another important variable and is maintained at the optimizer-specified value during the simulation by increasing inlet pressure. The lower bond is limited by the suction pressure of the cracked-gas compressor after the quenching section, which should not fall below atmospheric pressure to avoid any oxygen in-leak. A high value of steam to hydrocarbon ratio (SR) is favorable to lower coke formation and to increase run length but at the expense of increased steam consumption and pressure drop in the reactor. The inlet temperature (Tin) is bounded between 500 and 700 °C, which can be achieved in practice. Realistic constraints should be included in any optimization study. In the present study, two constraints were included: the first one is the maximum reactor wall temperature and the second one is run length. The upper limit of tube outer temperature at any point along the reactor length was fixed at 1100 °C. The run length of the cracker was limited between 30 and 76 days. These are based on typical industrial practice. 5. Results and Discussion Optimization of an industrial LPG thermal cracker is studied considering different, conflicting objective functions. The results are presented and discussed for biobjective optimization first and then for triobjective optimization. 5.1. Biobjective Optimization. Three cases of biobjective optimization were studied. Results for these cases are presented and discussed in this section. Case I. Maximization of Ethylene and Propylene Production per Year. The two objectivesstotal ethylene and propylene production per yearswere calculated assuming a decoking time of 36 h (as per industrial practice) and 8500 h of operation time per year, as follows. PC2) ) FC2)toperation
(7)
PC3) ) FC3)toperation
(8)
Here, toperation is the operation time given by toperation ) 8500 - td
(9)
where td is the decoking time in a year and can be calculated as follows:
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
8500 × 36 (10) RL × 24 where the run length (RL) is in days. Figure 2a shows the Pareto-optimal set obtained by NSGAII-aJG after 100 generations with 50 chromosomes, for maximizing annual ethylene and propylene production simultaneously. The generated Pareto set is reasonably smooth and well distributed. If we go from any one point to any other, say, from A to C in Figure 2a, ethylene production increases while propylene production decreases. Thus the solutions in Figure 2a comprise a Pareto front. The optimal values of two decision variables: Fin and Tin were nearly at constant at their respective upper bound (namely, 14 t/(h furnace) and 700 °C respectively), because maximization of both ethylene and propylene production flow rate are the goals. The optimal values of the other decision variabless COT, COP, and SR in Figure 2b-dsshow that the Pareto is due to their conflicting effect on the objectives. Figure 2b shows that COT has a strong opposing effect on PC2) and PC3); high COT maximizes ethylene production while low COT maximizes propylene production. For example, in the beginning of the Pareto front, low temperature (Figure 2b) and high pressure (Figure 2c) favor propylene production and so more propylene
9527
is produced. Ethylene production is increased by decreasing COP (partly due to increasing SR, Figure 2c and d) and increasing COT, in the later part of the Pareto front. Optimal SR is constant around 0.2 at ethylene production up to 32 kt/y but is scattered at higher values of ethylene production (Figure 2d). Table 3 compares the values objectives (PC2), PC3)), decision variables (Fin, COT, COP, SR, Tin), heat duty, and run length (RL) for three selected chromosomes A, B, and C in Figure 2a. It can be seen that total production of ethylene increases with COT but both propylene production and run length decrease. Furthermore, total heat duty increases with increasing COT, which is important because of cost of fuel and cracking coil. If the cracking coil works at high temperature, it has to be changed more frequently. Another key factor is run length. A high value of COT increases coke formation inside the coil, and the wall temperature increases rapidly requiring shutdown of the unit for decoking. These indicate that run length and total heat duty can be considered as additional objectives for the optimization of a thermal cracking process. These studies are presented in section 5.2. Case II. Maximization of Ethylene Selectivity Index and Minimization of Severity Index. The Pareto-optimal front
td )
Figure 2. (a) Optimal Pareto front for maximization of ethylene and propylene production per year. (b-d) Optimal values of decision variables corresponding to the Pareto front in part a. Table 3. Operating Conditions and Objective Values Corresponding to Chromosomes A, B, and C from Biobjective Optimization Case I shown in Figure 2 and Chromosomes G, H, and I from Triobjective Optimization Case B shown in Figure 8 selected chromosomes from biobjective optimization case I
selected chromosomes from triobjective optimization case B
quantity
A
B
C
G
H
I
Fin (t/(h furnace)) COT (°C) COP (bar) SR (kg steam/kg feed) Tin (°C) PC2) (kt/y) PC3) (kt/y) duty (Gcal/y) run length (d)
13.999 800.42 2.378 0.2068 698.56 30.375 22.794 624.970 58
13.998 820.76 1.812 0.2120 698.56 34.594 20.531 692.433 44
13.998 838.29 1.663 0.2464 697.59 37.439 18.389 744.221 30
13.965 801.60 2.334 0.3268 699.13 30.618 22.354 637.189 66
13.961 816.52 2.295 0.2433 699.13 33.734 20.929 677.414 40
13.953 838.95 1.844 0.393 697.56 37.600 18.135 761.165 30
9528
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
Figure 3. Results for simultaneous optimization of selectivity and severity indices (case II). (a) Optimal Pareto set. (b-d) Optimal values of decision variables corresponding to the Pareto front in part a. (e) Variation of constraint (run length) with selectivity index.
(Figure 3a) for the maximization of selectivity index and minimization of severity index is smooth and wide. Trends of the decision variables and constraint corresponding to this Pareto front (Figures 3b-e) are somewhat different from those in case I. For example, inlet flow rate and COP have high values in the beginning in order to minimize severity (more propylene production) whereas COT has a low value initially and steadily increases to increase the selectivity index. To increase the selectivity, the optimizer mainly increases the COT until its upper bound (Figure 3c) beyond which the feed flow rate is decreased (Figure 3b). When the COT increases, coke formation is boosted and consequently the COP and furnace run length decrease (Figure 3d-e). A low COP and high COT are favorable for β-scission reaction and hence for selectivity. The results in Figure 3 are consistent with this. The steam ratio and inlet temperature (not shown in Figure 3 for brevity) were at their upper and lower bound, respectively. High values of SR were selected, partly to maximize run length; note that high SR decreases coke formation and therefore the operation time is increased. A low inlet temperature is favorable for selectivity, run length, and propylene production because the average process gas temperature and coke formation will then be low.
Figure 4. Pareto-optimal set for maximization of ethylene production per year and minimization of heat duty per year (case III). The red squares are the calculated values of heat duty for the optimal Pareto set in case I, for comparison.
Case III. Maximization of Ethylene Production and Minimization of Heat Duty per Year. Figure 4 shows the Pareto set for maximization of annual ethylene production and minimization of annual heat duty. Optimal values of decision variables corresponding to this Pareto front are shown in Figure 5. As expected, heat duty increases with ethylene production
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
9529
Figure 5. Decision variables corresponding to the optimal Pareto set in Figure 4 for case III: maximization of ethylene production per year and minimization of heat duty per year. Table 4. Operating Conditions and Objective Values Corresponding to Chromosomes D, E, and F Shown in Figure 4 (Case III) quantity Fin (t/(h furnace)) COT (°C) COP (bar) SR (kg steam/kg feed) Tin (°C) PC2) (kt/y) PC3) (kt/y) duty (Gcal/y) run length (d)
chromosome D chromosome E chromosome F 10.554 841.00 1.555 0.2010 699.78 29.534 13.089 574.905 32
13.208 826.38 2.324 0.2040 699.96 33.535 18.593 661.468 30
13.874 840.46 1.503 0.2548 699.76 37.501 17.954 742.461 30
(Figure 4). Compared to the values on this Pareto front, the calculated annual heat duty for case I is higher as shown in Figure 4. This is not surprising because the optimizer is required to minimize heat duty. For the same reason, inlet temperature was chosen at its upper bound by the optimizer (not shown in Figure 5 for brevity). Table 4 shows the values of objectives, decision variables, and run length of three selected chromosomes from Figure 4 (case III). These values are different compared to those in Table 3. For example, the feed flow rate varies from its lower bound to its upper bound in case III (Table 4 and Figure 5a) whereas it is at the upper bound in case I (Table 3). Values of COT for chromosomes D, E, and F are higher than those for chromosomes A and B. It is clear from Tables 3 and 4 that heat duty is reduced in case III compared to case I by sacrificing propylene production. Also, run length in the case of chromosomes D, E, and F is shorter than that in case of chromosomes A and B. Therefore, triobjective optimization of ethylene and propylene production along with heat duty and/or run length is desirable. Two triobjective optimization cases are discussed in the next section. The trend of optimal values of decision variables is different from the previous cases. The feed flow rate is at its lower bound in the first half of the Pareto and then increases steadily to reach
the upper bound at high values of ethylene production (Figure 5a). This is because the optimizer increases COT first to increase ethylene production (Figure 5b). When the COT is near its upper bound, the feed flow rate is increased to increase ethylene production, which also results in increased heat duty. In this region, COT has a decreasing trend to achieve high selectivity but later increases, when the feed flow rate is near its upper bound, in order to increase conversion and hence annual ethylene production. The COP is high in the beginning (Figure 5c) because of low SR and high COT and then decreases with increasing of feed flow rate and COT. The steam ratio (Figure 5d) is practically constant at its lower bound because, for a particular COT, a high steam ratio needs more heat input. 5.2. Triobjective Optimization. Two cases of triobjective optimization are described in this section. Case A. Maximization of Ethylene Selectivity Index and Run Length, and Minimization of Severity Index. As discussed in case II in section 5.1, when optimization of two indices is considered, the run length of the reactor varied between its low and upper bounds (Figure 3e). Hence, triobjective optimization taking all three objectives was carried out. Figure 6 presents the Pareto-optimal solutions generated by maximizing ethylene selectivity (C2)/C2 ratio), run length, and minimization of severity (C3-/C3) ratio) simultaneously. It can be seen that when run length increases both severity and selectivity decrease. Hence, run length and selectivity are conflicting, whereas run length maximization and severity minimization do not generate nondominated solutions. It is interesting that the range of two objectives: selectivity and severity (Figure 6b) is the same as that in Figure 3a, Compared to biobjective results, the Pareto set from triobjective optimization has some scatter. Optimal values of the third objective function (run length) from triobjective optimization are compared in Figure 6c with the constraint (run length) values in case II. It can be seen that the former values of run length
9530
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
Figure 6. Pareto-optimal set obtained from the simultaneous optimization of ethylene selectivity, severity index, and run length: (a) three-dimensional plot and (b and c) plots of one objective against another. Values of objectives and run length in case II: (b and c) biobjective optimization of selectivity and severity indices, for comparison.
are better than those from biobjective optimization (case II). Hence, it can be concluded that triobjective optimization is necessary to improve the run length while optimizing selectivity and severity indices. Figure 7 shows the optimal values of decision variables corresponding to the Pareto front in Figure 6. Inlet temperature is at its lower bound and hence not shown in Figure 7. All decision variables have trends similar to those in case II (Figures 3b-d) except for some scatter, particularly in feed flow rate (Figure 7a). Hence, the overall discussion of trends in Figure 7 is similar to case II and not repeated here. Case B. Maximization of Ethylene and Propylene Production and Minimization of Heat Duty. Maximization of total ethylene production per year and minimization of total heat duty per year was studied earlier as case III. It was found that the total heat duty can be added as a third objective function to the two objectives in case I. Furthermore, minimization of heat duty reduces the fuel used and emission of gases (e.g., COx, SOx, and NOx) to the environment. The Pareto-optimal results obtained by triobjective optimization are plotted in two plots: total propylene production per year versus total ethylene production per year (Figure 8a) and total heat duty per year versus total ethylene production per year (Figure 8b). Figure 8a shows that annual ethylene and propylene production from triobjective optimization covers a wide range and appears scattered in the two-dimensional plots. Several points in the top-right corner are nearly comparable to the Pareto of case I: biobjective optimization of ethylene and propylene. The heat duty for the Pareto set in case I was calculated and shown in Figure 8b. Similarly, the annual propylene production for case III was calculated and included in Figure 8a. It can be see
from these plots that the calculated values for the Pareto set in case I are comparable to the optimal values of heat duty in case B. However, the Pareto set for case B has a wider range. In other words, the optimizer generated more diverse solutions for triobjective optimization. The Pareto set in case III is superior to that in case B and calculated values of heat duty in case I (Figure 8b). However, in the first part of this plot, the Pareto set for cases B and III are nearly same. Further, it can be seen from Figure 8a that the calculated values of annual propylene production for case III in this region are less than the optimal values found in case B. The dependency of the Pareto-optimal set in Figure 8 on the decision variables (Figure 9) can be explained in a similar fashion as for the results of biobjective optimization (Figures 4 and 5). There is more scatter in Figure 9 compared to Figure 5, and change in the feed flow rate and COT at PC2) in the range 28 to 30 kt/y is very sharp in Figure 9. Optimal values of objectives and operating conditions for selected chromosomes from biobjective optimization case I are compared with those for chromosomes from triobjective optimization case B in Table 3. Chromosomes G, H, and I in case B are selected such that their annual ethylene and propylene production are comparable to that for chromosomes A, B, and C in case I, respectively. However, the total heat duty is different in each case, which may be partly due to differences in annual ethylene and propylene production rates. 6. Conclusions MOO of an industrial LPG thermal cracker was carried out using NSGA-II-aJG and a detailed, first principles model.
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
9531
Figure 7. Decision variables corresponding to the Pareto set in Figure 6 for the triobjective optimization of ethylene selectivity, severity, and run length.
range of 30-38 and 17.5-22.5 kt/y for ethylene and propylene respectively. The best operating point from the Pareto-optimal solutions can be determined by considering user’s preferences, demand for products, process operability, safety analyses, etc. The variation of values of decision variables with the objectives at the optimum in all cases can be explained qualitatively, which shows that MOO results obtained by NSGA-II-aJG are reliable. Coil outlet temperature is the most important decision variable in all cases studied, whereas the influence of other decision variables such as feed flow rate, COP, SR, and inlet temperature depends on the objectives considered. Results of this study show that triobjective optimization of the LPG thermal cracker is desirable to find a wider range of optimal solutions compared to those by biobjective optimization, for the objectives considered. Appendix A: Reactor Model Equations The governing equations for a steam cracking reactor are as follows.29 The mass balance for each of the molecular and free radical species is given by
( )(∑ )
dFj πdin2 ) dz 4
Figure 8. Pareto-optimal set for simultaneous maximization of annual propylene and ethylene production, and minimization of heat duty per year (case B): (a) Total propylene production per year versus total ethylene production per year and (b) total heat duty versus total ethylene production per year. Results of biobjective optimization (case I, ethylene and propylene production and, case III, ethylene production and heat duty) including calculated values of total propylene production for Pareto set in case III and total heat duty for Pareto set in case I are shown for comparison.
Several MOO problems with two and three objectives were solved successfully by NSGA-II-aJG. The Pareto set for optimizing annual production of ethylene and propylene has the
Rijri
(A-1)
i
Here, Fj is the flow rate of j ) 1, 2, 3,..., the total number of molecular and free radical species (excluding steam), the summation on the right-hand side is for all reactions, and the rate of the ith reaction, ri, is expressed as
( )∏
ri ) ko,i exp where Cj )
(
-Ei RT
Fj (Fsteam
Cj
)
P ( ) RT +∑F) j
j
(A-2)
j
(A-3)
9532
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009
Figure 9. Decision variables corresponding to the Pareto front in Figure 8 for the maximization of annual ethylene and propylene production and minimization of total heat duty per year.
Note that the product term on the right-hand side of eq A-2 and Cj are only for the reactants (both molecular species and free radicals) taking part in the ith reaction. See the Nomenclature for the definition of all other symbols. The energy balance is written considering heat of reactions and heat flux:
Coke formation at any time during the thermal cracking of feed and at any position in the reactor is calculated from the continuity equation:21
(∑ F C
n The total rate of coke formation is expressed as Rc ) ∑i)1 rc,i where n is the number of coking reactions producing coke. The number of coke precursors and coking reactions contributing to coke formation can be found in the work of Towfighi et al.22
)(dTdz ) ) Q(z)πd
j p,j + FsteamCp,steam
j
out -
( )(∑ πdin2 4
)
ri(-∆Hi) (A-4)
i
T ∫298
∆Cp,i dT. The first summation is with ∆Hi ) ∆Hi,298 + over all molecular species whereas the summation on the righthand side is for all cracking reactions. In this study, the heat flux profile, Q(z), is initialized based on the independent simulation of the fire-box and later manipulated during the reactor simulation to achieve the specified COT. The momentum balance equation is formulated to accommodate the additional pressure drop in the bends by using a friction factor, Fr.
(
)
( )
Pt dPt d 1 1 1 1 dT ) + Fr (A-5) + 2 MmPt RG RT dz dz Mm Mm T dz
(
)
The friction factor for straight tubes and bends is respectively given by Re-0.2 din
(A-6)
Re-0.2 ζ + din πRb
(A-7)
Fr ) 0.092 Fr ) 0.092 where Re )
)(
din dinG Λ 0.051 + 0.19 and ζ ) 0.7 + 0.35 µm 90° Rb
(
)
∂Cc(t, z) ) Rc ∂t
(A-8)
Nomenclature Cj ) concentration of species j in gas mixture, kmol/m3 Cp,j ) specific heat capacity of species j, kJ/(kmol K) Cc ) coke concentration C2)/C2 ) selectivity index, ethylene to ethane fraction in the product stream, % wt/% wt C3-/C3) ) severity index, sum of propane, propylene, propadiene, ethane, ethylene, acetylene, methane, and hydrogen fraction to propylene fraction, % wt/% wt COP ) coil outlet pressure, bar COPst ) set point of coil outlet pressure, bar COT ) coil outlet temperature, °C COTst ) set point of coil outlet temperature, °C din ) tube internal diameter, m dout ) tube external diameter, m Ei ) activation energy of reaction i, kcal/kmol Fin ) feed flow rate, T/(h furnace) Fj ) molar flow rate of species j other than steam, kmol/s Fsteam ) molar flow rate of steam, kmol/s FT ) total molar flow rate of gas mixture including steam, kmol/s Gm ) mass flux of the gas mixture, kg/(m2 s) ko,i ) frequency factor for reaction i, 1/s or m3/(kmol s) L ) total length of reactor, m laJG ) length of the replacing jumping gene
Ind. Eng. Chem. Res., Vol. 48, No. 21, 2009 lchrom ) total length (number of binaries) of a chromosome lsubstr ) length of substring Mm ) molecular weight of gas mixture, kg/kmol Npop ) number of chromosomes in the population Ngen ) number of generations p ) pressure, bar Pin ) inlet pressure, bar pc ) crossover probability PC2) ) total ethylene production per year, kt/y PC3) ) total propylene production per year, kt/y PQ ) total heat duty per year, Gcal/y pm ) mutation probability paJG ) jumping probability for the JG operator Q(z) ) heat flux at length z of the reactor, kcal/(m2 s) R ) ideal gas constant, 8.31451 J/(mol K) Rb ) bend radius, m Re ) Reynolds number ri ) rate of reaction i, kmol/(m3 s) rc,i ) reaction rate of coke formation, kg/(m3 s) RL ) run length, days Sr ) seed for random number generator SR ) steam (to hydrocarbon) ratio, kg steam/kg feed t ) time, h tinc ) time increment for calculation of coke thickness, day T ) temperature, K Tin ) inlet temperature, °C Tw ) reactor wall temperature, °C Tw,max ) maximum reactor wall temperature, °C z ) distance from the reactor entrance, m Greek Symbols ∆Ηi ) heat of reaction i, kJ/kmol ∆z ) reactor length increment, m Λ ) angle described by the bend, taken to be 180° R ) unit conversion factor Rij ) stoichiometric coefficient of species j in reaction i µm ) viscosity of gas mixture, kg/(m s)
Literature Cited (1) Tagliabue, M.; Ferrari, M.; Pollesel, P. Increasing value from steam cracker olefin streams. Pet. Technol. Q. 2004, 9, 145–149. (2) Grantom, R. L.; Royer, D. J. Ethylene. In Ullmann’s Encyclopedia of Industrial Chemistry, 5th ed.; Wiley-VCH: Weinheim , 1997. (3) Eliceche, A. M.; Patrica Hoch, N. C.; Brignole, E. A. Optimal operation of an ethylene plant at variable feed conditions. Comput. Chem. Eng. 1995, 19, S223. (4) Edwin, E. H.; Balchen, J. G. Dynamic optimization and production planning of thermal cracking operation. Chem. Eng. Sci. 2001, 56, 989. (5) Towfighi, J.; Karimaei, M.; Karimzadeh, R. Application of Dynamic Programming Method in optimal control of light hydrocarbons pyrolysis reactors. Sci. Iran. 1998, 5 (3), 148. (6) Shahrokhi, M.; Nedjati, A. Optimal temperature control of propane thermal cracking reactor. Ind. Eng. Chem. Res. 2002, 41, 6572. (7) Schulz, E. P.; Bandoni, J. A.; Diaz, M. S. Optimal Shutdown Policy for Maintenance of Cracking Furnaces in Ethylene Plants. Ind. Eng. Chem. Res. 2006, 45, 2748. (8) Bhaskar, V.; Gupta, S. K.; Ray, A. K. Application of multiobjective optimization in chemical engineering. ReV. Chem. Eng. 2000, 16 (1), 1. (9) Nandasana, A. D.; Ray, A. K.; Gupta, S. K. Applications of the non- dominated sorting genetic algorithm (NSGA) in chemical reaction engineering. Int. J. Chem. Reaction Eng. 2003, 1, R21 (available at http:// www.bepress.com/ijcre/). .
9533
(10) Masuduzzaman; Rangaiah, G. P. Multi-Objective Optimization Applications in Chemical Engineering. In Multi-ObjectiVe Optimization: Techniques and Applications in Chemical Engineering; Rangaiah, G. P., Ed.; World Scientific: Singapore, 2009. (11) Tarafder, A.; Lee, B. C. S.; Ray, A. K.; Rangaiah, G. P. Multiobjective optimization of an Industrial ethylene reactor using a nondominated sorting genetic algorithm. Ind. Eng. Chem. Res. 2005, 44, 124. (12) Li, C.; Zhu, Q.; Geng, Z. Multiobjective particle swarm optimization hybrid algorithm: an application on industrial cracking furnace. Ind. Eng. Chem. Res. 2007, 46, 3602. (13) Gao, X.; Chen, B.; He, X.; Qiu, T.; Li, J.; Wang, C.; Zhang, L. Multi-objective optimization for the periodic operation of the naphtha pyrolysis process using a new parallel hybrid algorithm combining NSGAII with SQP. Comput. Chem. Eng. 2008, 32, 2801. (14) Guria, C.; Bhattacharya, P. K.; Gupta, S. K. Multiobjective optimization of reverse osmosis desalination units using different adaptations of non-dominated sorting genetic algorithm (NSGA). Comput. Chem. Eng. 2005, 29, 1977. (15) Froment, G. F.; Van de Steen, B. O.; Van Damm, P. S.; Narayanan, S.; Goossens, A. G. Thermal cracking of ethane and ethane-Propane mixtures. Ind. Eng. Chem. Process Des. DeV. 1976, 15, 495. (16) Sundaram, K. M.; Froment, G. F. Modeling of thermal cracking Kinetics. 1. Thermal cracking of ethane, propane and their mixtures. Chem. Eng. Sci. 1977, 32, 601. (17) Sundaram, K. M.; Froment, G. F. Modeling of thermal cracking Kinetics. 3. Radical mechanisms for the pyrolysis of simple paraffins, olefins and their mixtures. Ind. Eng. Chem. Fundam. 1978, 17, 174. (18) Herriott, G. E.; Eckert, R. E.; Albright, L. Kinetic of propane pyrolysis. AIChE J. 1972, 18, 84. (19) Van Damm, P. S.; Narayanan, S.; Froment, G. F. Thermal cracking of propane and propane-propylene mixtures: Pilot plant versus industrial data. AIChE J. 1975, 21, 1065. (20) Sundaram, K. M.; Froment, G. F. kinetics of coke deposition in the thermal cracking of propane. Chem. Eng. Sci. 1979, 34, 635. (21) Plehiers, P. M.; Reyniers, G. C.; Froment, F. G. Simulation of the run length of an ethane cracking furnace. Ind. Eng. Chem. Res. 1990, 29, 636. (22) Towfighi, J.; Niaei, A.; Karimzadeh, R.; Saedi, G. Systematic and modeling representations of LPG thermal cracking for olefin production. Kor. J. Chem. Eng. 2006, 23 (1), 8. (23) Deb, K. Multi-objectiVe Optimization Using EVolutionary Algorithms; Wiley: Chichester, UK, 2001. (24) Deb, K.; Pratap, A.; Agarwal, S.; Meyarivan, T. A. Fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans. EVol. Comput. 2002, 6, 181. (25) Kasat, R. B.; Gupta, S. K. Multiobjective optimization of an industrial fluidized-bed catalytic cracking unit (FCCU) units using genetic algorithm (GA) with the jumping genes operator. Comput. Chem. Eng. 2003, 27, 1785. (26) Agrawal, N.; Rangaiah, G. P.; Ray, A. K.; Gupta, S. K. Multiobjective optimization of the operation of an Industrial low density polyethylene tubular reactor using genetic algorithm and its jumping gene adaptations. Ind. Eng. Chem. Res. 2006, 45, 3182. (27) Agrawal, N.; Rangaiah, G. P.; Ray, A. K.; Gupta, S. K. Design stage optimization of an Industrial low density polyethylene tubular reactor for multiple objectives using NSGA-II and its jumping gene adaptations. Chem. Eng. Sci. 2007, 62, 2346. (28) Van Geem, K. M.; Reyniers, M. F.; Marin, G. B. Two Severity Indices for Scale-Up of Steam Cracking Coils. Ind. Eng. Chem. Res. 2005, 44, 3402. (29) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design; Wiley & Sons: New York, 1990.
ReceiVed for reView September 19, 2008 ReVised manuscript receiVed January 16, 2009 Accepted January 22, 2009 IE801409M