Multiobjective Optimization of a Fixed Bed Maleic Anhydride Reactor

Article pubs.acs.org/IECR

Multiobjective Optimization of a Fixed Bed Maleic Anhydride Reactor Using an Improved Biomimetic Adaptation of NSGA-II Pranava Chaudhari and Santosh K. Gupta* Department of Chemical Engineering, Indian Institute of Technology, Kanpur 208016, India ABSTRACT: A fixed bed maleic anhydride reactor is simulated using an improved reactor model. A simplified Langmuir−Hinshelwoodtype reaction kinetics is used, and the kinetic as well as the reactor parameters (for a full-scale pilot plant reactor and an industrial reactor) are tuned. The model is then used to solve a few optimization problems involving a single and multiple objectives from among maximum productivity, minimum operating cost, and minimum pollution. This is carried out using NSGA-II-aJG. A semiempirical procedure is suggested to reduce the scatter in the Pareto optimal solutions obtained for the three-objective optimization problems, and replace solutions which are associated with extreme sensitivity (this is different from obtaining a robust Pareto set, where all the optimal points have reduced sensitivity). Thus, our focus is twofold: modeling and multiobjective optimization of MA reactors, and improving the algorithm to deal with the extreme sensitivity of the three-objective problem (which may be of benefit in other problems as well). In addition, it is observed that a more recent algorithm, Alt-NSGA-II-aJG, biomimicking the altruism of honeybees, converges to the optimal solutions faster than does NSGA-II-aJG for a two-objective optimization problem, but it gives inferior solutions for a three-objective problem.

Table 1. Details of the Catalyst and the Reactors (PPR and IndR)a

1. INTRODUCTION Maleic anhydride (MA) is an extremely versatile chemical. It is an excellent cross-linking agent because of three active sites (two carboxyl groups and a double bond). Commercial production of MA started in the 1930s with benzene as a raw material in a multitube fixed bed reactor. After the mid-1980s, most of the industries switched over to n-butane (Bu) feed due to environmental pressure on the use of benzene and lower Bu prices. The first commercial plant using this feed was Monsanto’s J. F. Queeny plant in Illinois, started in 1974. MA is used mainly in the production of unsaturated polyester resins. Other applications include agricultural chemicals, maleic and fumaric acids, lube oil additives, and a number of other specialty chemicals and organic intermediates. The global production of MA in 2010 was approximately 1.7 million metric tons. Its consumption is expected to grow an average of 5.6% per year from 2010 to 2015 and 3.5% per year from 2015 to 2020. The conversion of Bu to MA involves complex kinetics on the internal surface of the porous catalyst. This has been studied by several investigators.1−5 Bej and Rao6,7 proposed an eight-parameter redox kinetic model which they tuned using their data on an isothermal differential reactor. Sharma et al.8 have used experimental data9 on a nonisothermal full-scale pilot plant reactor (PPR) using a VPO catalyst and the reactor model proposed by Wellauer et al.,9 to tune a simplified Langmuir− Hinshelwood-type six-parameter kinetic model. In contrast to developments in the reaction and reactor engineering of MA production, studies on the optimization of fixed bed MA reactors are scarce. Early research on the optimization of these reactors focused on process modifications,10 either by improving the catalyst or by using alternative reactors, e.g., fluidized bed reactor,11,12 circulating fluidized bed reactor,13 membrane reactor,14 dual-coolant (salt-bath) system,9 and twozone fluidized bed reactor.15 However, fixed bed reactors are still in operation in several industries and a few studies on their optimization have been reported.9,16 The yield of MA was maximized in these studies. Unfortunately, a higher yield (eq A48 in Appendix 1) © 2012 American Chemical Society

reactor and catalyst properties 0

T (K) D (m) H (m) ρcat (kg/(m3 of porous catalyst)) εcat rcat (m) ξ (Å) τ E1 (kJ kmol−1) E2 (kJ kmol−1) E3 (kJ kmol−1) Rg′ (kJ kmol−1 K−1) Rg (m3 atm kmol−1 K−1) K2 (atm−1) α1 α3 p n yO0 2 0 0 0 yH0 2 , yMA , yCO , yCO 2 0 yN2 a

values

ref

453.15 0.025 5 1400 (PPR), 1600 (IndR) 0.402 0.0015 637 3 93 100 155 000 93 100 8.31 0.082 06 310 0.54 0.54 1 5.5 0.21 0.0 0 1 − yO0 2 − yBu

8, 9 8, 9 8, 9 8 24 8, 9 24 8 8 8 8

8 8 8 8 8 8

Values of the IndR are the same as for the PPR, unless otherwise specified.

0 , and the of MA is associated with a lower inlet Bu flow rate, FBu maximization of the yield leads to low overall production rates, FMA, of MA. This suggests the optimization of these reactors using multiple objectives (multiobjective optimization, MOO), to give a better picture for decision making.

Received: Revised: Accepted: Published: 3279

October 4, 2011 January 10, 2012 January 25, 2012 January 25, 2012 dx.doi.org/10.1021/ie202276q | Ind. Eng. Chem. Res. 2012, 51, 3279−3294

Industrial & Engineering Chemistry Research

Article

In eqs 1b and 1c, p and n are stoichiometric coefficients depending on the porous catalyst used. The values used in the present study are 1 and 5.5 for p and n, respectively, as proposed by Sharma et al.8 The simple one-dimensional heterogeneous model19 proposed by Sharma et al.8 and modified slightly is used. This accounts for the exothermic fast reactions. Radial and axial dispersion for both heat and mass transfer are assumed to be negligible. There are a total of seven components (1, Bu; 2, O2; 3, H2O; 4, MA; 5, CO2; 6, CO; 7, N2) present in the system. Mole balance equations are written for each of the seven components for both the gas and the catalyst phases. In contrast to the present detailed model, Sharma et al.8 have written mole balances for only a few components and use stoichiometry and the rates of reaction to solve for the remaining. They have used FT, the molar flow rate of the gas, as a constant. It has been reported by Sharma et al.8 that the temperature within the catalyst is almost uniform and varies with axial location only. We have also assumed this to be so since the Prater number for the catalyst is quite low.20 The overall heat transfer coefficient from the coolant (salt) to the packed catalyst bed is obtained using the correlation given by Beek,21 with an empirical multiplying factor,22 f Uov. This differs from Sharma et al.8 Appendix 1 gives the complete set of coupled equations used in the present study (a full discussion of the model details are not being provided here for the sake of brevity, since they can be inferred from Sharma et al.8). Some of the equations (for the gas phase) are ordinary differential equations of the initial value kind (ODE-IVP), while some (for the catalyst phase) are ODEs of the boundary value kind (ODE-BVP). We have applied orthogonal collocation23 (with four collocation points using the matrices for the symmetrical case for a sphere) to convert the ODE-BVPs into a set of nonlinear algebraic equations. The cylindrical catalyst pellet actually used has been treated as an equivalent sphere (eq A4). The system of resulting ODE-IVPs is stiff23 because of the highly reactive region near the inlet of the reactor. The ODE-IVPs (using the Gear method,23 code D02EJF, in the NAG library) and the system of nonlinear algebraic equations (using the modified Powell hybrid algorithm, code C05NBF in the NAG library) are solved simultaneously. 2.2.1. Tuning of Parameters. Pilot Plant Reactor (PPR). Experimental data on a PPR have been provided8,9 at five sets of operating conditions, referred to as PPR cases 1−5. Three of these (cases 1, 3, and 5) are used for fitting some of the parameters of the present model. The rate constants, kj (j = 1, 2, 3), of the three reactions in eqs 1a−1c are written in terms of their values, kj673 (j = 1, 2, 3), at a reference temperature of 673 K (eq A15). The latter are “tuned” (best-fit values obtained) by fitting the experimental data for cases 1, 3, and 5 of the PPR. Two additional reactor parameters, namely, the bed porosity,

In this study, we first improve the reactor model of Sharma et al.8 and then tune and validate the model parameters using experimental data.8,9 Several single objective optimization (SOO) and two- and three-objective problems are then solved. The performance of two recent biomimicked adaptations17,18 of the elitist MOO algorithm, NSGA-II, are also compared for one of the twoobjective optimization problems solved, as well as one of the threeobjective optimization problems solved to see the performances of these algorithms.

2. REACTOR MODELING 2.1. Reaction Scheme and Mathematical Model. The reaction mechanism involved in MA production by the partial oxidation of Bu comprises of both consecutive and parallel reactions. The simplified reaction scheme proposed by Sharma et al.8 and used here is given by k1

C4 H10 (Bu) + 3.5O2 → C4H2O3 (MA) + 4H2O (desired reaction)

(1a)

k2

C4 H2O3 + pO2 → (6 − 2p)CO + (2p − 2)CO2 + H2O

(decomposition reaction)

(1b)

k3

C4 H10 + nO2 → (13 − 2n)CO + (2n − 9)CO2 + 5H2O

(total oxidation reaction)

(1c)

Table 2. Operating Conditions for the Different Cases of PPR8 and the IndR9 operating conditions 0

−2 −1

G (kg m s ) 0 yBu (mol %) PT0 (atm) TS (K)

case 1 (PPR)

case 2 (PPR)

case 3 (PPR)

case 4 (PPR)

case 5 (PPR)

IndR

1.24 1.82 1.9 673.15

1.24 1.86 1.65 663.15

1.19 0.75 1.64 653.15

1.19 0.75 1.64 643.15

1.21 1.26 1.62 623.15

1.3 1.8 2.5 648.15

Table 3. Best-Fit (Tuned) Values of the Parameters parameter k1673 (kmol/(kg of porous catalyst·s·atm0.54)) k2673 (kmol/(kg of porous catalyst·s)) k3673 (kmol/(kg of porous catalyst·s·atm0.54)) εb (m3 of gas mixture/m3 of bed) f Uov

tuned values (PPR)

tuned values (IndR)

Sharma et al.8

7.5 × 10−6

7.5 × 10−6

0.96 × 10−6

2.5 × 10−5

2.5 × 10−5

0.29 × 10−5

8.4 × 10−7

8.4 × 10−7

1.5 × 10−7

0.42

0.40

−

1.1

0.71

−

Table 4. Best Values of the Computational Parameters

a

problem

populn size

max generation

0 lstring (G0, yBu , PT0 , T0)

lstringa (TS/Zopt; TS,1, TS,2)

pcross

pmute

pjump

ljump

problems 1, 2 (PPR and IndR) problems 3, 4 (PPR) (NSGA-II-aJG) problem 4b (PPR) (Alt-NSGA-II-aJG) problems 3, 4 (IndR) problem 5 problem 6 (NSGA-II-aJG) problem 6c (Alt-NSGA-II-aJG)

40 60 100 100 60 60 60

80 300 80 300 400 600 600

− 15 20 20 20 10 15

10 30 30 30 30 10, 10, 10 30, 30, 30

0.3 0.95 0.90 0.95 0.95 0.90 0.95

2/10 3/90 3/90 3/110 3/110 3/70 3/150

0.3 0.4 0.4 0.4 0.4 0.2 0.2

5 5 15 5 10 5 5

lstring for Zopt, TS,1, and TS,2 are only for problem 6. bNqueen = 8. cNqueen = 6. 3280

dx.doi.org/10.1021/ie202276q | Ind. Eng. Chem. Res. 2012, 51, 3279−3294


Article

εb (not reported by Sharma et al.8), and f Uov, used in eq A38 with the Beek21 estimate of the overall heat transfer coefficient, are also tuned. Results for the other two cases, cases 2 and 4, are used for validation of the tuned model. The details of the reactor and the catalyst used, as well as the operating conditions used, are given in Tables 1 and 2.8,9,24 Tuning of the five parameters is done by minimizing the cumulative weighted sum-of-squares error between the modelpredicted results and experimental values. NSGA-II-aJG17 is used for this single objective optimization (SOO) problem, with the objective function given by

changed from the value of 1400 kg/m3 porous catalyst used for the PPRs (see Table 1) to 1600 kg/m3 porous catalyst for the IndR. This was done to obtain a better fit of the simulated results,9 while still being in the acceptable range8 of values. A procedure similar to that used for the PPR is used to tune the two parameters. The only difference is that the weighting factor, wR, for the CO/CO2 ratio was taken as zero since simulation results for the CO/CO2 ratio were not reported9 for the IndR.

3. OPTIMIZATION 3.1. Safety Aspects of the Reactor. The maximum temperature inside the catalytic packed bed reactor should be below about 713.15−723.15 K to prevent catalyst deactivation.25 A maximum value of 713.15 K is used in the present study. In addition, the temperature of the gas mixture in the reactor system should lie below the autoignition temperature (AIT) outside the packed bed (flame propagation is minimal inside the bed since the catalyst packing acts as a flame arrester).8 Figure 1 shows that the AIT varies with the fraction of Bu in the air as well as with the pressure.26 An additional tem-

Min I(k1673 , k2673 , k3673 , ε b, fU ) ov ⎡⎧ NT (j) ⎫ NCase ⎢ p ⎛ TModel(i , j) − TExp(i , j) ⎞2 ⎪ ⎪ ⎟⎟ ⎬ ≡ ∑ ⎢⎨ ∑ wT ⎜⎜ TModel(i , j) ⎪ ⎝ ⎠ ⎪ ⎢ j = 1 ⎣⎩ i = 1 ⎭ ⎛ XBu,Model(j) − XBu,Exp(j) ⎞2 ⎟ + wX ⎜⎜ ⎟ XBu,Exp(j) ⎝ ⎠ ⎛ SMA,Model(j) − SMA,Exp(j) ⎞2 ⎟ + wS⎜⎜ ⎟ SMA,Exp(j) ⎝ ⎠ ⎛ RModel(j) − RExp(j) ⎞2 ⎟ + wR ⎜⎜ ⎟ RExp(j) ⎝ ⎠ ⎤ ⎛ ΔP Model(j) − ΔP Exp(j) ⎞2 ⎥ ⎟ ⎥ + wΔP ⎜⎜ ⎟ ΔP Exp(j) ⎝ ⎠ ⎥ ⎦

(2)

In eq 2, NCase (=3) is the total number of cases used for fitting the data (PPR cases 1, 3 and 5), NTp(j) is the total number of data points used to fit the temperature at different locations in case j, T(i,j) is the experimental (subscript “Exp”) or modelpredicted (subscript “Model”) value of the temperature at location i for case j, XBu(j) is the experimental or modelpredicted value of the conversion of n-butane (eq A46) for case j, SMA(j) is the value of the selectivity (eq A47) for MA production for case j, R(j) is the value of the CO/CO2 ratio at the exit of the reactor for case j, and ΔP(j) is the value of the pressure drop through the reactor bed for case j. The weighting factors, wT (=40), wX (=2), wS (=2), wR (=1), and wΔP (=1), have been used to provide proper weightings to the different errors. 2.2.2. Industrial Reactor (IndR). Wellauer et al.9 provide simulation results (and not actual data points) for an IndR. The details of the reactor and the catalyst used, as well as the operating conditions used, are given in Tables 1 and 2. This reactor has also been tuned in the present study. We have taken the values of the three kinetic parameters (kj673; j = 1, 2, 3) to be the same as in the case of the PPR (Table 3), but we have tuned only the two other (unreported) reactor parameters, εb and f Uov, for the IndR. Values of the latter parameters may be different in the IndR compared to those in the PPR because the bed porosity depends on the packing strategy, which may differ for the IndR, and the heat transfer coefficient in a multitubular reactor may be smaller than that estimated for a single tube used in the case of the PPR (the same correlation of Beek21 is used for both cases in this study). Another parameter, ρcat, was

Figure 1. Experimental values (thin lines) of AIT vs fraction of butane in air. The maximum allowable temperature, Tcnstrnt, used in the present study also shown by a bold line.

perature constraint, T < Tcnstrnt, is used both at the inlet and the exit of the packed bed. Tcnstrnt is a function of the Bu concentration and is taken empirically (see Figure 1) as if (100FBu/FT) < 1.3;

Tcnstrnt = 740 − 80(FBu/FT)

if (100FBu/FT) ≥ 1.3;

Tcnstrnt = 520

(3)

Figure 1 shows the Tcnstrnt used in this study, about 50 K below the value for 3 atm, the upper bound of pressure in the present study. 3.2. Formulation of Optimization Problems. An important objective in the optimization of a multitubular industrial MA reactor is to maximize the total (exit) flow rate of MA from the reactor at minimum cost. The latter involves both the fixed cost and the operating cost. The exit flow rate of MA from the reactor can be increased either by increasing the 0 number of tubes at the same input Bu flow rate, FBu , per tube 0 (increases the fixed cost), or by increasing FBu while keeping the number of tubes constant (increases the operating cost), or both. Only the latter is an option in an operating reactor. The 0 0 value of FBu depends on G0, yBu , and PT0 , the superficial mass velocity, the mole fraction of Bu in the feed, and the inlet 3281



Article

Figure 2. Comparison (cases 1, 3, 5) and validation (cases 2, 4) of model-predicted results with experimental8 (Exp) values on the PPR; (a) temperature profiles, (b) Bu conversion, (c) selectivity, (d) CO/CO2 ratio at the end of the reactor, and (e) pressure drop. “Our Sim” represents the predictions of the tuned model from the present study, while “Sharma et al. Sim” represents the simulation results of Sharma et al.8

pressure, respectively. There are two other parameters, T0 and TS, the feed temperature and the salt-bath coolant temperature, respectively, which determine the temperature inside the reactor tube and also play an important role in the optimal operation of the reactor. These five can be taken as decision variables in an optimization study. However, one must solve a few simpler SOO and MOO problems first to develop intuition and help obtain solutions for more complex MOO problems,27 as done in several earlier real-life industrial optimization studies.27 We start with an SOO problem for an operating reactor with a 0 , using the inlet values G0 = 1.24 kg m−2 s−1, fixed value of FBu 0 = 0.0186, yO0 2 = 0.21, and PT0= 1.65 atm (PPR, case 2) and yBu other values given in Table 1. Two decision variables, TS and T0, are used first, but it is observed that the results are quite insensitive to the latter. Hence, only TS is used as the decision variable (with T0 = 453.15 K; PPR, case 2). The SOO problem is described by

problem 1 (SOO): Max: I(TS) ≡ FMA (kmol/s)

(4a)

subject to (s.t.)

bounds: 573.15 ≤ TS ≤ 713.15 K

(4b)

constraints: T ≤ Tcnstrnt

(4c)

where Tcnstrnt at inlet and exit:

inside the reactor:

eq 3

Tcnstrnt = 715 K

(4d) (4e)

NSGA-II-aJG17 is used for solving this problem. 3282



Article

problem 2 (MOO):

A simple two-objective optimization problem similar to eqs 4a−4e 0 (G0 = 1.24 kg m−2 s−1, yBu = 0.0186, yO0 2 = 0.21, PT0 = 1.65 atm, and T0 = 453.15 K; PPR, case 2) is then solved, which maximizes the MA production while simultaneously minimizing the undesired waste products, CO and CO2:

Max: I1(TS) ≡ FMA (kmol/s)

(5a)

Min: I2(TS) ≡ FCO + FCO2 (kmol/s)

(5b)

subject to (s.t.) bounds and constraints: eqs 4b−4e

(5c)

Problems 1 and 2 are also solved for the industrial reac0 tor (IndR), with G0 = 1.3 kg m−2 s−1, yBu = 0.018, yO0 2 = 0.21, T0 = 453.15 K, and PT0 = 2.5 atm, and using the other values in Table 1. Having built up some skills in obtaining optimal solutions, more complex and relevant SOO and MOO problems are 0 solved. FBu (operating cost) is now allowed to be changed 0 through the use of four more parameters, G0, yBu , PT0 , and T0, as decision variables along with TS. The SOO problem maximizes 0 FMA, as in problem 1. In the two-objective problem, FBu is minimized while FMA is maximized simultaneously. This reduces the operating cost while maximizing the profit. The two problems are described as problem 3 (SOO): 0 , P 0 , T 0, T ) ≡ F Max: I(G0 , yBu T S MA (kmol/s)

(6a)

bounds:

Figure 3. Comparison of model-predicted results from the tuned model; ―, present study; ---, Wellauer et al.,9 for the IndR.

⎛ kg ⎞ 1.0 ≤ G0 ≤ 4.0 ⎜ 2 ⎟ ⎝ m ·s ⎠

(6b)

0 ≤ 0.04 (mol %) 0.01 ≤ yBu

(6c)

1.5 ≤ P T0 ≤ 3.0 (atm)

(6d)

Figure 4. Solution of the optimization problem 2 (MOO), for both (a, b) the PPR case 2 and (c, d) the IndR. Result for problem 1 (SOO) also shown by ◆ for both. 3283



Article

Figure 5. (a) Optimal Pareto set for the MOO problem 4 for the PPR: ▽, using NSGA-II-aJG; △, using Alt-NSGA-II-aJG. ■, Optimal solution for SOO problem 3; ◆, experimental values8 for the five cases. (b−f) Decision variables corresponding to different points of the Pareto set for MOO, optimal point for SOO, and simulation points of the five cases for the PPR.8,9

453.15 ≤ T 0 ≤ 673.15 (K)

(6e)

573.15 ≤ TS ≤ 713.15 (K)

(6f)

problem 5 (MOO):

s.t. constraints: eqs 4c−4e

(8a)

0 0 Min: I2(G0 , yBu , P T0 , T 0 , TS) ≡ FBu (kmol/s)

(8b)

(6g) 0 Min: I3(G0 , yBu , P T0 , T 0 , TS) ≡ FCO + FCO2 (kmol/s)

problem 4 (MOO):

(8c)

0 Max: I1(G0 , yBu , P T0 , T 0 , TS) ≡ FMA (kmol/s)

(7a)

0 0 Min: I2(G0 , yBu , P T0 , T 0 , TS) ≡ FBu (kmol/s)

(7b)

s.t. bounds and constraints:

eqs 6b−6g

(8d)

The best values of the computational parameters used for all cases are given in Table 4. The reactions involved in this process are highly exothermic. This leads to high temperatures near the inlet and to poor performance of the reactor as well as to a possibility of catalyst deactivation.9 Wellauer et al.9 have suggested the use of a dual salt-bath system (dividing the reactor cooling system into two

s.t. bounds and constraints:

eqs 6b−6g

0 Max: I1(G0 , yBu , P T0 , T 0 , TS) ≡ FMA (kmol/s)

(7c)

The final optimization problem solved involves three objective functions, with the additional minimization of the exit flow rates of CO + CO2: 3284



Article

Figure 6. (a) Optimal Pareto set for the MOO problem 4 for the IndR. ■, Optimal solution for the SOO problem 3; ◆, reported9 simulation point for MOO. (b−f) Decision variables corresponding to different points of the Pareto set for MOO, optimal point for SOO, and simulation point for the IndR.

segments with different coolant temperatures, TS,1 and TS,2, with the switchover occurring at axial location, Zopt) to provide improved performance in the reactor. A three-objective optimization problem is solved for this case involving additional decision variables:

573.15 ≤ TS,1 ≤ 713.15 (K)

(9e)

573.15 ≤ TS,2 ≤ 713.15 (K)

(9f)

0.5 ≤ Zopt ≤ 3.0 (m)

(9g)

problem 6 (MOO): 4. RESULTS AND DISCUSSION 4.1. Tuning of Parameters. Experimental data are available for the temperature at a few locations in the PPR as well as for the conversion of Bu (eq A46), the selectivity (eq A47), the CO/CO2 ratio at the exit, and the total pressure drop for cases 1−5. These are shown as data points (“Exp”) in Figure 2. Data on three cases (1, 3, and 5) are used to tune the five model parameters, k1673, k2673, k3673, εb, and f Uov. The best-fit values of the parameters for the PPR so obtained are given in Table 3. We obtain higher values for the kinetic parameters than reported.8 This is due to the modifications in the model used in the present study. Figure 2 gives a comparison of the model-predicted results with the experimental values for all five cases, three of which have been used for tuning (comparison) as well as the two cases which have not been so used, but used for validation. Figure 2 also gives the simulation results8,9 for

0 Max: I2(G0 , yBu , P T0 , T 0 , TS,1, TS,2 , Zopt)

≡ FMA (kmol/s)

(9a)

0 Min: I1(G0 , yBu , P T0 , T 0 , TS,1, TS,2 , Zopt) 0 ≡ FBu (kmol/s)

(9b)

0 Min: I3(G0 , yBu , P T0 , T 0 , TS,1, TS,2 , Zopt)

≡ FCO + FCO2 (kmol/s)

(9c)

s.t. bounds and constraints: eqs 6b−6e

(9d) 3285



Article

Figure 7. (a) I1 (in increasing order) for the PPR (problem 5); corresponding values of (b) I2 and (c) I3; (d−h) corresponding values of the decision variables. ▼, Original (optimal) solutions; ▲, improved solutions.

results for the IndR obtained using the fitted parameters and the simulation results9 is shown in Figure 3. Reasonable agreement is observed. It may be added that our attempts to obtain better agreement for the temperature did not succeed (we did not change the heats of reactions from their values for the PPR). 4.2. Optimization. The solutions of problem 2 (MOO), for both the PPR case 2 and the IndR, are shown in Figure 4. An increase in FMA is accompanied by an increase in FCO + FCO2, indicating that the curve in Figure 4 is a Pareto set. The corresponding values of the salt-bath temperature, TS, for each of the Pareto points are also given in Figure 4. In order to get a higher FMA, the single decision variable, TS, should be increased. However,

the conversion of Bu, selectivity, CO/CO2 ratio at the exit, and overall pressure drop. It is observed that our tuned model matches the experimental data slightly better than the earlier model.8 Interestingly, the factor, f Uov, used for the heat transfer correlation is close to unity, indicating that Beek’s21 correlation used for the single-tube PPR is appropriate. The tuned values of the parameters for the IndR are also given in Table 3. We observe that the bed porosity, εb, for the IndR is not too different from that for the PPR, but the multiplying factor, f Uov, for the heat transfer coefficient is smaller than for the PPR, as expected. The comparison of the simulation 3286



Article

Figure 8. (a) I1 (in increasing order) for the PPR (problem 6; using NSGA-II-aJG); corresponding values of (b) I2 and (c) I3.

Figure 9. (a) I1 (in increasing order) for the PPR (problem 6) using (△) Alt-NSGA-II-aJG and (▼) NSGA-II-aJG; corresponding values of (b) I2 and (c) I3.

since the value of E2, the activation energy for the decomposition reaction (eq 1b) is larger than E1, the value for the desired reaction (eq 1a), this leads to a simultaneous increase in FCO + FCO2. Figure 4 also shows the solution of problem 1 (SOO). The highest value of FMA is indicated corresponding to the constraint eq 4e on the temperature inside the reactor (higher than TS). The optimal Pareto set for problem 4 (MOO) along with the optimal point for problem 3 (SOO) are shown in Figures 5a and 6a for the PPR and the IndR, respectively. The length of binaries used for TS is about twice that used for the other

decision variables in the chromosome (see Table 4). The use of equal number of binaries for each of the decision variables leads to poorer results. A few points from these Pareto sets were chosen randomly and the value of TS was varied manually, with the values of the other decision variables fixed. It was observed that the points in Figures 5 and 6 were, indeed, the best. This confirms the correctness of the MOO code used. The computational parameters were also varied around their best values given in Table 4 to study their effects. Increasing the population size worsened the optimal values of the important decision 3287



Article

variable, TS (while not affecting the Pareto set), whereas a decrease in the population size led to a reduction in the range of the Pareto set. Making pjump = 0 (i.e., using NSGA-II without the JG adaptation) again led to a similar worsening of TS. An increase in the length, ljump, of the jumping gene, to 10 leads to an increase in the scatter of the results for all the decision variables. It is observed that the solutions of the SOO problem 3 (SOO) lie near the extreme right corner of the optimal Pareto sets for the MOO problem 4, in the cases of both the PPR and the IndR. The experimental values8 for cases 1−5 of the PPR are also shown in Figure 5. It is observed that cases 1−4 lie close to the Pareto set, while case 5 does not (operating the pilot plant under near-optimal conditions was not planned by the authors). The operation of the IndR9 is, similarly, not optimal (Figure 6 a). The optimal values of the five decision variables associated with the different points on the Pareto set are also shown in

Figures 5b−f and 6b−f. An interesting phenomenon is observed, associated with the complex constraints in eqs 4d and 4e. The first is associated with the conditions at the inlet and exit of the reactor, while the second is related to the maximum temperature permissible inside the reactor. Below FMAof about 2.5 × 10−7 kmol/s, for both the PPR and the IndR, 0 yBu increases to a limiting value of 0.013 a few times and then drops as FMA increases. This is dictated by eq 4d and Figure 1. Simultaneously, the optimal value of TS increases until it attains its maximum value (corresponding to eq 4e) at FMA of about 2.5 × 10−7 kmol/s. The values of G0 and PT0 exhibit a sudden 0 jump every time yBu drops in the initial phase so as to compensate for its fall, and to achieve the next point of the Pareto set. Above FMA of about 2.5 × 10−7 kmol/s, TS decreases as FMAincreases, while the reactor temperature always lies at its 0 maximum permissible value and yBu increases continuously. The maximum achievable value of TS is higher in the case of the PPR than in the case of the IndR because heat transfer in the latter case (f Uov = 0.71) is lower than in the former (f Uov = 1.1). For the three-objective problem 5, the chromosomes in the final generation are rearranged in increasing order of I1 (≡FMA). The optimal solutions, I1−I3, are plotted as a function of the rearranged chromosome number. Arranging solutions in this manner28 has certain advantages, one of which is that the scatter can be observed easily (scatter is more difficult to identify in three-dimensional plots of the optimal objective functions). Figure 7a−c indicates that these solutions, indeed, comprise a Pareto set because I1 improves (increases) while the other two objective functions, I2 and I3, worsen (decrease). It is observed that the optimal value of I1 in problem 5 corresponding to the same I2 is lower than that for the twoobjective problem 4. The best set of computational parameters for this problem is given in Table 4. Increasing the population size leads to more scatter in the optimal results of I2 and I3, while decreasing the population size leads to a reduction in the range of the Pareto set. Increasing pjump leads to more scatter in the plots of I2 and I3, while decreasing pjump leads to a reduction in the range of the Pareto set. An decrease in ljump leads to an increase in the scatter of I2 and I3. An interesting observation for the three-objective problem is that there are four points (shown by filled inverted triangles in Figure 7b) which do not lie close to the neighboring points (particularly for I2), even with the best values of the computational parameters. The ε-constraint method29 is used to confirm that these four points are, indeed, optimal. A detailed study reveals that there is an acute sensitivity associated with these four optimal points, with I2 worsening signif icantly with a slight improvement of I1 and I3. This is not so for the other optimal points. A similar phenomenon of acute sensitivity and

Figure 10. Flow chart for Alt-NSGA-II-aJG.18

Table 5. Information on Specific Heat, Standard Enthalpy, and Diffusion Volumesa component i b

Bu O2 H2O MAb CO2 CO N2

αCp,i (kJ kmol−1 K−1)

βCp,i (kJ kmol−1 K−2)

106γCp,i (kJ kmol−1 K−3)

107δCp,i (kJ kmol−1 K−4)

10−5Hf,i298.15 (kJ kmol−1)

diffusion vol (∑m νm)i

55.16 25.74 30.38 −13.05 26.02 26.88 27.31

0.182 0.0130 0.00962 0.348 0.0435 0.0069 0.0052

0.0 −3.860 1.180 218.07 −14.840 −0.820 −0.0042

0.0 0.0 0.0 4.832 0.0 0.0 0.0

−1.265 0.0 −2.42 −3.983 −3.935 −1.105 0.0

86.7 16.3 13.1 68.25 26.9 18 18.5

Unless otherwise indicated, αCp,i, βCp,i, γCp,i, and Hf,i298.15 are from Becker24 and (∑m νm)i are from Poling et al.34 bαCp,i, βCp,i, γCp,i, δCp,i, and Hf,i298.15 values are from Perry and Green.35 a

3288



Article

scatter of the optimal Pareto solutions was observed in our earlier three- and four-objective optimization problems in chemical engineering,17,28,30−32 but was not followed up. Such sensitivity of optimal solutions is not acceptable in real-life systems, and these points need to be replaced by other, less sensitive optimal points. The values of TS for these four chromosomes are found to be quite low compared to those for the neighboring points (see Figure 7) and are artif icially increased so that they come closer to their neighbors. The shifted points so obtained (filled triangles) are also optimal, as confirmed using the ε-constraint technique. The corresponding optimal 0 solutions show less scatter (including values of yBu ) and sensitivity. The remaining scatter in the plots of the decision variables can be reduced further using a similar procedure, but this is not done. This semiempirical procedure of “doctoring” results obtained from NSGA-II-aJG is recommended for computationally intense real-life MOO problems exhibiting sensitivity, at least until some better methodology is developed (since there are plenty of optimal points, with no large gaps, in the Pareto sets, it may not matter much in the present case. However, there may be other cases where this procedure may be of use). The optimal values of the several decision variables, associated with the Pareto solutions, and shown in Figure 7d−h, show trends similar to those for problem 4 (except for TS being lower). The procedure suggested for the replacement of some solutions which are associated with extreme sensitivity is slightly different from that suggested by Ramteke and Gupta33 for generating a robust Pareto, where all the optimal points have reduced sensitivity (are robust). The solutions of the three-objective problem 6 for the PPR, with two salt-bath temperatures, TS,1 and TS,1, and a switchover at Zopt, are shown in Figure 8a−c, arranged in the same manner as Figure 7a−c. Interestingly, better results are obtained when an equal number of binaries is used for all the decision variables (see Table 4). The computational parameters were varied around their best values to study their effects. Increasing or decreasing the population size led to more scatter in the results. Making pjump = 0 (i.e., using NSGA-II without the JG adaptation) led to a reduction in the range of the Pareto set, while increasing pjump to 0.4 led to more scatter in the results. An increase in pcross to 0.95 led to more scatter in the results, while an increase in pmute to 3/lchr led to more scatter. An increase in the length, ljump, of the jumping gene to 10 led to a reduction in the range of the Pareto set as well as to an increase in the scatter of the results. Ten chromosomes were found to lie far from the neighboring points in the Pareto set. The semiempirical approach described for problem 5 is used to improve the results. Figure 8 shows the results after this doctoring up. It would be expected27 that the scatter in the Pareto front increases as the number of decision variables or objective functions increases. The trends exhibited by the optimal values of the decision variables for problem 6 (results not shown for brevity, but can be supplied on request) are similar to those in problems 4 and 5. We observe from Figure 8 that, for the first few chromosomes, I1 and I2 improve while I3 worsens. In contrast, for other chromosomes, I1 improves while I2 and I3 worsen. Plotting of the optimal results for problems involving three or more objectives in the manner done in Figures 7 and 8 helps in observing this kind of behavior. A comparison of the results of problems 5 and 6 shows that the use of a dual salt-bath system leads to higher values of FMA 0 and FCO + FCO2 for the same FBu . This means that the dual salt-

bath system is superior in terms of production of MA, but at the cost of higher pollution. Wellauer et al.9 had suggested that the dual salt-bath system is superior since the production of MA is higher, but since he did not consider the third objective, he failed to realize that it is associated with a cost. The usefulness of solving the three-objective problems is evident. The two-objective optimization problem 4 and the threeobjective optimization problem 6 are also solved using the more recent algorithm,18 Alt-NSGA-II-aJG, biomimicking the altruistic behavior of honeybees. The computational parameters for Alt-NSGA-II-aJG have been varied over a wide range: for problem 4, pj = 0.0−0.5, pcross = 0.9−0.95, pmute = 1/lchr−4/lchr, Nqueen = 6−10; for problem 6, pj = 0.0−0.3, pcross = 0.9−0.95, pmute = 1/lchr−3/lchr, Nqueen = 4−8, and the best values identified. These are given in Table 4. The Pareto optimal solutions are shown in Figure 5 (problem 4) and Figure 9 (problem 6). For problem 4, it is observed that Alt-NSGA-IIaJG gives the same (superposable) Pareto solutions in almost half the number (8000) of total function calculations than those taken by NSGA-II-aJG (18 000). For the three-objective problem 6, the same number (36 000) of function calculations has been used for Alt-NSGA-II-aJG as for NSGA-II-aJG (for which the results are shown in Figure 8). It is observed (Figure 9) that Alt-NSGA-II-aJG gives a smaller range of solutions and is also associated with a higher amount of scatter than the results from NSGA-II-aJG (Figure 8). In fact, our attempts to obtain faster convergence or better results for the three-objective problem 6 using Alt-NSGA-II-aJG (by varying the computational parameters over very wide ranges) failed. Interestingly, Ramteke and Gupta18 also observed that Alt-NSGA-II-aJG performed better than NSGA-II-aJG for the two-objective optimization of a phthalic anhydride reactor (they did not, however, study a three-objective optimization problem for the phthalic anhydride reactor using either Alt-NSGA-II-aJG or NSGA-II-aJG). It is too early at this stage to generalize this observation (superior performance of Alt-NSGA-II-aJG over NSGA-II-aJG for the two-objective problem and its inferior performance for a three-objective problem).

4. CONCLUSIONS A fixed bed maleic anhydride reactor is simulated using an improved reactor model. The model parameters are tuned using some results in the literature for pilot plant reactors and an industrial reactor. The model is used to solve some optimization problems involving a single and multiple objectives from among maximizing productivity, minimizing operating cost, and minimizing pollution. Two recent, computationally efficient biomimicked adaptations, NSGA-II-aJG and AltNSGA-II-aJG, are used to generate results. It is found that the latter performs better than NSGA-II-aJG for a two-objective problem, but does not perform as well for a three-objective problem. A semiempirical procedure is suggested for the replacement of some solutions for the three-objective optimization problem, which are associated with extreme sensitivity. This also reduces the scatter in the Pareto optimal solutions. This new procedure is different from that suggested by Ramteke and Gupta33 for generating a robust Pareto, where all the optimal points have reduced sensitivity (are robust). The resulting improved solutions would be of more use for industrial applications. 3289



■

Article

⎛ ⎞ kmol * λ = k3p α3 ⎜ R3, Bu, λ ⎝ kg of porous cat ·s ⎟⎠

APPENDIX 1: COMPLETE SET OF MODEL EQUATIONS

Overall Mole Balance

(A14)

⎡⎛ ⎤ ⎞ ⎢⎜ Ei ⎟⎜⎛1 − 673 ⎟⎞⎥ ; k j = k 673 exp j ⎢⎣⎜⎝ 673R g′ ⎟⎠⎝ T ⎠⎥⎦

dFi = −A(1 − ε b)ρcat9 i; dZ

j = 1, 2, 3

i = 1, Bu; 2, O2 ; 3, H2O; 4, MA; 5, CO2 ; 6, CO; 7, N2

(A15)

pi , λ = Ci ,cat|λ R gT

(A1)

1

∫λ= 0 {λ2R̅ i , λ} dλ

9i = 3

Energy Balance (A2)

Mole Balance inside the Catalyst Particle

′ 2R̅ i , λρcat rcat 1 d ⎡ 2 dCi ,cat|λ ⎤ 2 λ ≡ ∇ | = C ; ⎢ ⎥ spherical i ,cat λ dλ ⎦ Die λ2 dλ ⎣ i = 1, 2, ..., 7

at λ = 0:

dλ dCi ,cat|λ= 1

at λ = 1: ≡−

dλ

=0

(A3b)

λ≡

where

⎧ 3 ⎫ ⎪ ⎪ 2 * ⎨ ∑ [( −ΔHj , λ)R j , λ]λ ⎬ dλ Q=3 λ= 0 ⎪ ⎪ ⎩ j=1 ⎭

(A18)

(A3d)

Nc

2

FT =

since lcat = rcat

(A20)

∑ Fi i=1

(A4)

Effective Diffusivity

*λ RBu, ̅ λ = R1,* λ + R3,

(A5)

* λ + nR3, *λ R̅ O2, λ = 3.5R1,* λ + pR2,

ε Die = cat Di ; τ

(A6)

* λ − 5R3, *λ R̅ H2O, λ = −4R1,* λ − R2,

(A7)

* λ + R2, *λ RMA, ̅ λ = − R1,

(A8)

* λ − (2n − 9)R3, *λ R̅ CO2, λ = −(2p − 2)R2,

(A9)

* λ − (13 − 2n)R3, *λ R̅ CO, λ = −(6 − 2p)R2,

(A10)

R̅ N2, λ = 0

(A11)

* λ = k2 R2,

⎛ ⎞ kmol ⎜ ⎟ (1 + K2pMA, λ )2 ⎝ kg of porous cat ·s ⎠

i = 1, 2, ..., 7

(A21)

(For eq A21, see eq 5.183, Elnashaie and Elshishini.20) 1 1 1 = + ; Di D k, i Dm, i

D k, i = (97 × 10−10)ξ

i = 1, 2, ..., 7 (A22)

T ; Mw, i

i = 1, 2, ..., 7 (A23)

7 yj 1 1 ; = ∑ Dm, i 1 − yi j = 1 Di , j

α1 pBu, λ

⎛ ⎞ kmol R1,* λ = k1 ⎜ ⎟ 1 + K2pMA, λ ⎝ kg of porous cat ·s ⎠

(A19)

′ = 2rcat d pe = 2rcat

spherical equivalence: ′ rcat rcat lcat r = = cat ; 2 3 3 2rcatlcat + 2rcat

1

(G0)2 (1 − ε b) dP T 1 =− dZ ρ0d pe ε b3 1.01325 × 105 ⎡ 150(1 − ε )μ ⎤ 0 b gas ⎥ P T T FT ×⎢ + 1.75 0 ⎢ ⎥ PT T 0 F 0 G d pe T ⎣ ⎦

(A3c)

r ′ rcat

(A17)

Ergun Equation

= ∇spherical Ci ,cat|λ= 1

′ ⎛ k g, ircat F ⎞ ⎜Ci ,cat|λ= 1 − P T i ⎟ Die ⎜⎝ R gT FT ⎟⎠

Q (A(1 − ε b)ρcat) − Uov πD(T − TS) dT = N dZ ∑i =comp 1 FiCp , i

∫

(A3a)

dCi,cat|λ= 0

(A16)

i = 1, 2, ..., 7

j≠i

(A12)

pMA, λ

yi = (A13) 3290

Ci ,cat|λ ; Nc ∑i = 1 Ci ,cat|λ

(A24a)

i = 1, 2, ..., 7 (A24b)



Article

⎛ m2 ⎞ ⎟⎟ Di , j ⎜⎜ ⎝ s ⎠ =

Heats of Reaction (All H’s at Temperature T)

ΔH1, λ = HMA, λ + 4HH2O, λ − HBu, λ − 3.5HO2, λ (1.43 × 10−7)T1.75

(A33)

(1.01325P T)Mi , j1/2[(∑m νm)i1/3 + (∑m νm)j1/3 ]2

i , j = 1, 2, ..., 7

; ΔH2, λ = (6 − 2p)HCO, λ + (2p − 2)HCO2, λ + HH2O, λ − HMA, λ − pHO2, λ

(A25)

(A34)

34

(For eq A25, see Poling et al. )

ΔH3, λ = (13 − 2n)HCO, λ + (2n − 9)HCO2, λ

−1 ⎡⎛ ⎞⎤ ⎛ 1 ⎞ ⎜ 1 ⎟⎥ ⎢ ⎟+ ; Mi , j = 2 ⎜⎜ ⎢⎣⎝ Mw, i ⎟⎠ ⎜⎝ Mw, j ⎟⎠⎥⎦

+ 5HH2O, λ − HBu, λ − nHO2, λ

i , j = 1, 2, ..., 7 (A26a)

Molar Specific Heats and Enthalpies of Formation of Different Components

(∑ vm)i = sum of the atomic diffusion volumes of the mth

T ⎛ kJ ⎞ ⎟; Hi , λ(T ) = Hf,298.15 + CT p , i dT ⎜⎝ i kmol ⎠ T = 298.15

∫

m

atoms of component i

(A26b)

i = 1, 2, ..., 7

Gas−Film Mass Transfer Coefficient, kg,i (Dwivedi and Upadhyay36)

Shi =

Shi =

k g, id pe film Dm, i

Regas =

Sci ,gas =

ρ0 =

i = 1, 2, ..., 7

(Values of αCp,i, βCp,i, δCp,i, and Hf,i298.15are given in Table 5.) Calculation for Overall Heat Transfer Coefficient (Beek;21 Hyman22)

i = 1, 2, ..., 7 (A28)

Uov = fU Uov,Beek ov ⎧ Kg ⎪ = fU ⎨ 0.4 [2.58Regas0.33Prgas0.33 ⎪ ov d pe ⎩

Gd pe μgas

(A29)

μgas film ρgasDm, i

P T0 R gT 0

(A37)

(A27)

;

P T0 ρgas = ρ0 T0 PT T

(A36)

2 3 ⎛ kJ ⎞⎟ CT ; p , i = αC p, i + βCp , iT + γCp , iT + δCp , iT ⎜⎝ kmol ·K ⎠

0.4548 (Regas)0.5931(Sci ,gas)1/3 ; εb

i = 1, 2, ..., 7

(A35)

;

⎫ ⎪ ⎛ kJ ⎞ + 0.094Regas0.8Prgas0.4]⎬ ⎟ ⎜ ⎪ 2 ⎭ ⎝ m ·s·K ⎠

i = 1, 2, ..., 7 (A30)

FT0 FT

⎛ kJ kg CT p ,mix μgas ⎜ kg·K m·s ⎜ Prgas = kJ Kg ⎜⎜ m s·K · ⎝

(A31a)

0 Mw,mix

⎞ ⎟ ⎟ ⎟⎟ ⎠

(A38)

(A39)

(A31b)

CT p ,mix =

Nc

F0 0 Mw,mix = ∑ i0 Mw, i F i=1 T

Nc

∑ i=1

(A31c)

T Fi Cp , i FT Mw,mix

Nc

Viscosity of Air; Sutherland’s Formula (Crane37)

Mw,mix =

i=1

3/2

T + 120 ⎛ T ⎞ μgas = μref ref ⎜ ⎟ T + 120 ⎝ Tref ⎠

∑

(Pa·s)

(A32b)

μref = 1.827 × 10−5 Pa·s

(A32c)

⎛ kg ⎞ Fi ⎟ Mw, i ⎜ ⎝ kmol ⎠ FT

(A40)

(A41)

K g = 0.024035(1 + 0.00317(T − 273.15)

(A32a)

Tref = 291.15 K

⎛ kJ ⎞ ⎜ ⎟ ⎝ kg·K ⎠

⎛ J ⎞ ⎟ − 0.0000021(T − 273.15)2 ) ⎜ ⎝ m·s·K ⎠

(A42)

(For eq A42, see Kannuluik and Carman.38) 3291



Article T Cp,i = molar heat transfer capacity of component i at temperature T, kJ/(kmol·K) T Cp,mix = specific heat transfer capacity of gas mixture at temperature T, kJ/(kg·K) dpe = equivalent catalyst particle diameter, m D = diameter of reactor, m Die = effective diffusivity of component i in the gas mixture at (Z, λ), m2/s Dk,i = Knudsen diffusivity of component i in gas mixture at (Z, λ), m2/s Di,j = molecular diffusivity of component i in component j Dm,i = diffusivity of component i in the gas mixture at (Z, λ), m2/s film Dm,i = diffusivity of component i in the gas film around the catalyst particle at Z Ej = activation energy for jth reaction, kJ/kmol; j = 1, 2, 3 f Uov = multiplying factor for Uov Fi = molar flow rate of component i at Z, kmol of i/s FT = total molar flow rate of gas mixture at Z, kmol of i/s G = superficial mass velocity, kg of gas/(m2 superficial·s) H = height of the reactor, m Hi,λ = enthalpy of component i at (Z, λ), kJ/kmol ΔHj,λ = enthalpy of jth reaction at (Z, λ), kJ/kmol kj = rate constant for jth reaction, kmol/(kg of porous catalyst·s·atm−αj) kj 673 = rate constant at 673 K for jth reaction, kmol/(kg of porous catalyst·s·atm−αj) kg,i = mass transfer coefficient of component i in the gas film around catalyst particle at Z, m/s Kg = thermal conductivity of air at temperature T, kJ/ (m·s·K) lchr = length of binaries for one (full) chromosome lstring = length of binaries for single variable ljump = length of jumping gene Mw,i = molecular weight of component i, kg/kmol Mw,mix = molecular weight of the gas mixture at (Z, λ), kg/kmol Nc = number of components in gas mixture (=7) Nqueen = number of queen bees in a population pcross = probability of crossover pmute = probability of mutation pjump = probability of jumping gene pi,λ = partial pressure of component i at (Z, λ), atm PT = pressure at height, Z, atm r = radial location inside spherical catalyst particle, m rcat = radius of the cylindrical catalyst particle, m rcat ′ = equivalent radius of the spherical particle, m 9 i = overall rate of consumption of component i, per kilogram of porous catalyst at Z, kmol of i/kg of porous catalyst R̅ i, λ = rate of consumption of component i, per kilogram of porous catalyst basis at (Z, λ), kmol of i/kg of porous catalyst Rj,λ * = rate of jth reaction at (Z, λ), kmol of i/kg of porous catalyst Rg = gas constant, atm·m3/kmol·K Rg′ = gas constant, kJ/(kmol·K) Reg = Reynolds number at Z Sci,gas = Schmidt number for component i, at Z Shi = Sherwood number for component i, at Z SMA = overall selectivity for the desired product, MA, % T = temperature of bed at height Z, K TS = salt (coolant) temperature, K

Initial Conditions for ODEs

At Z = 0 Fi = Fi0 = yi0 FT0;

i = 1, 2, ..., 7

T = T0

(A44)

P T = P T0 FT0 =

(A43)

G0 0 Mw,mix

(A45a)

A (A45b)

Conversion, XBu, of Bu

XBu =

0 FBu − FBu ·100% 0 FBu

(A46)

Selectivity, SMA

4FMA ·100% 4FMA + FCO + FCO2 F = 0 MA ·100% FBu − FBu

SMA =

(A47)

Yield, YMA, (mol %)

■

F conversion· selectivity YMA = MA ·100 = 0 100 FBu

(A48)

APPENDIX 2: BRIEF DESCRIPTION OF TWO RECENT BIOMIMICKED ADAPTATIONS OF NSGA-II

NSGA-II-aJG

The following steps are incorporated while generating new daughter chromosomes after crossover and mutation.27,33 • Select a chromosome (sequentially) from the daughter population; check if the aJG operation is needed, using pjump. • If yes, select one location on this chromosome randomly between positions 1 and lchr − ljump. • The ljump binaries starting from the selected location are to be replaced by a new set of randomly generated binaries. Alt-NSGA-II-aJG.18

Figure 10 gives the Alt-NSGA-II-aJG algorithm.

■ ■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

NOMENCLATURE A = inside cross-sectional area of the reactor bed, m2 Ap = external surface area of catalyst particle, m2 Ci,cat|λ = moles of component, i, per cubic meter micropore volume, mol of i/m3 gas [not mol of i/m3 (pore + solid)] 3292



Article

(13) Contractor, R. M.; Bergna, H. E.; Horowitz, H. S.; Blackstone, C. M.; Chowdhry, U.; Sleight, A. W. Butane Oxidation to Maleic Anhydride in a Recirculating Solids Reactor. Catalysis 1988, 38, 645. (14) Mallada, R.; Menéndez, M.; Santamaría, J. Use of Membrane Reactors for the Oxidation of Butane to Maleic Anhydride under High Butane Concentrations. Catal. Today 2000, 56, 191. (15) Gascón, J.; Téllez, C.; Herguido, J.; Menéndez, M. Fluidized Bed Reactors with Two-Zones for Maleic Anhydride Production: Different Configurations and Effect of Scale. Ind. Eng. Chem. Res. 2005, 44, 8945. (16) Bej, S. K.; Rao, M. S. Selective Oxidation of n-Butane to Maleic Anhydride. 1. Optimization Studies. Ind. Eng. Chem. Res. 1991, 30, 1819. (17) Kasat, R. B.; Gupta, S. K. Multi-Objective Optimization of an Industrial Fluidized-Bed Catalytic Cracking Unit (FCCU) using Genetic Algorithm (GA) with the Jumping Genes Operator. Comput. Chem. Eng. 2003, 27, 1785. (18) Ramteke, M.; Gupta, S. K. Biomimicking Altruistic Behavior of Honey Bees in Multi-Objective Genetic Algorithm. Ind. Eng. Chem. Res. 2009, 48, 9671. (19) de Wasch, A. P.; Froment, G. F. Heat Transfer in Packed Beds. Chem. Eng. Sci. 1972, 27, 567. (20) Elnashaie, S. S. E. H., Elshishini, S. S. Modelling, Simulation and Optimization of Industrial Fixed Bed Catalytic Reactors; Gordon and Breach: London, 1993. (21) Beek, J. Advances in Chemical Engineering; Academic Press: New York, 1962; Vol. 3. (22) Hyman, M. Simulate Methane Reformer Reactions. Hydrocarbon Process. 1968, 47, 131. (23) Gupta, S. K. Numerical Methods for Engineers, 2nd ed.; New Age International: New Delhi, 2010. (24) Becker, C. Katalytische Wandreaktorkonzepte für MSA-Synthese und Methanol-Dampfreformierung. Ph.D. Thesis, Universität Stuttgart, Germany, 2002. (25) Alonso, M.; Lorences, M. J.; Pina, M. P.; Patience, G. S. Butane Partial Oxidation in an Externally Fluidized Bed-Membrane Reactor. Catal. Today 2001, 67, 151. (26) Chandraratna, M. R.; Griffiths, J. F. Pressure and Concentration Dependences of the Autoignition Temperature for Normal Butane + Air Mixtures in a Closed Vessel. Combust. Flame 1994, 99, 626. (27) Ramteke, M.; Gupta, S. K. Multi-Objective Genetic Algorithm and Simulated Annealing with Jumping Gene Adaptations. In MultiObjective Optimization: Techniques and Applications in Chemical Engineering; Rangaiah, G. P., Ed.; World Scientific: Singapore, 2008; p 91. (28) 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. (29) Haimes., Y. Y. Hierarchical Analysis of Water Resources Systems; McGraw-Hill: New York, 1977. (30) Guria, C.; Bhattacharya, P. K.; Gupta, S. K. Multi-Objective Optimization of Reverse Osmosis Desalination Units using Different Adaptations of the Non-Dominated Sorting Genetic Algorithm (NSGA). Comput. Chem. Eng. 2005, 29, 1977. (31) 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. (32) Khosla, D. K.; Gupta, S. K.; Saraf, D. N. Multi-Objective Optimization of Fuel Oil Blending using the Jumping Gene Adaptation of Genetic Algorithm. Fuel Process. Technol. 2007, 88, 51. (33) Ramteke, M.; Gupta, S. K. Biomimetic Adaptations of GA and SA for the Robust MO Optimization of an Industrial Nylon-6 Reactor. Mater. Manuf. Processes 2009, 24, 38. (34) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids, 5th ed.; McGraw-Hill: New York, 2001.

TS,i = salt (coolant) temperature of ith segment of dual saltbath system, K Uov = overall heat transfer coefficient, kJ/(m2·s·K) XBu = overall conversion of Bu, % yi = fraction of component i YMA = overall yield for the desired product, MA, % Z = axial location in the reactor, m Zopt = axial location of the reactor where the switchover is occurring for dual salt-bath system, m Superscript/Subscript

0 = feed values Greek Symbols

εb = void fraction in the catalyst bed (macrovoid), m3 of macrovoid/m3 of bed (porous catalyst + macrovoid) εcat = volume fraction of micropores in the catalyst pallet, m3 of micropores/m3 of porous catalyst λ = dimensionless radius of catalyst (=r/rcat ′ ) μgas = viscosity of the gas mixture at Z, Pa·s ξ = average pore size of catalyst, Å ρcat = density porous catalyst, kg of porous catalyst/m3 of porous catalyst ρgas = density of gas mixture in the macrovoid, kg of gas/m3 of gas ρgas = density of gas mixture at Z, kg of gas mixture/m3 of macropores τ = tortuosity of the catalyst pallet

■

REFERENCES

(1) Wohlfahrt, K.; Hofmann, H. Kinetik der Synthese von Maleinsäureanhydrid aus n-Butan. Chem. Ing. Tech. 1980, 52, 811. (2) Sharma, R. K.; Cresswell, D. L.; Newson, E. J. Kinetics and Reactor Modeling in Fixed-Bed Pilot Plant Production of Maleic Anhydride by Oxidation of n-Butane. Presented at the AIChE Meeting, San Francisco, 1984; Paper 69b. (3) Buchanan, J. S.; Sundaresan, S. Kinetics and Redox Properties of Vanadium Phosphate Catalysts for Butane Oxidation. Appl. Catal. 1987, 26, 211. (4) Centi, G.; Trifiro, F.; Ebner, J. R.; Franchetti, V. M. Mechanistic Aspects of Maleic Anhydride Synthesis from C4 Hydrocarbons over Phosphorus Vanadium Oxide. Chem. Rev. 1988, 88, 55. (5) Uihlein, K. Butanoxidation an VPO―Wirbelschichtkatalysatoren. Ph.D. Thesis, Fakultät für Chemieingenieurwesen, Universität Karlsruhe, Germany, 1993. (6) Bej, S. K.; Rao, M. S. Selective Oxidation of n-Butane to Maleic Anhydride. 2. Identification of Rate Expression for the Reaction. Ind. Eng. Chem. Res. 1991, 30, 1824. (7) Bej, S. K.; Rao, M. S. Selective Oxidation of n-Butane to Maleic Anhydride. 3. Modeling Studies. Ind. Eng. Chem. Res. 1991, 30, 1829. (8) Sharma, R. K.; Cresswell, D. L.; Newson, E. J. Kinetics and FixedBed Modeling of Butane Oxidation to Maleic Anhydride. AIChE J. 1991, 37, 39. (9) Wellauer, T. P.; Cresswell, D. L.; Newson, E. J. Optimal Policies in Maleic Anhydride Production through Detailed Reactor Modelling. Chem. Eng. Sci. 1986, 41, 765. (10) Ballarini, N.; Cavani, F.; Cortelli, C.; Ligi, S.; Pierelli, F.; Trifirò, F.; Fumagalli, C.; Mazzoni, G.; Monti, T. VPO Catalyst for n-Butane Oxidation to Maleic Anhydride: A Goal Achieved, or a Still Open Challenge? Top. Catal. 2006, 38, 147. (11) Contractor, R. M.; Sleight, A. W. Maleic Anhydride from C-4 Feedstocks using Fluidized Bed Reactors. Catal. Today 1987, 1, 587. (12) Stefani, G.; Budi, F.; Fumagalli, C.; Suciu, G. D. Fluidized Bed Oxidation of n-Butane: A New Commercial Process for Maleic Anhydride. Stud. Surf. Sci. Catal. 1990, 55, 537. 3293



Article

(35) Perry’s Chemical Engineering Handbook, 8th ed.; Perry, R. H., Green, D. W., Eds.; McGraw-Hill: New York, 2008. (36) Dwivedi, P. N.; Upadhyay, S. N. Particle-Fluid Mass Transfer in Fixed and Fluidized Beds. Ind. Eng. Chem. Process Des. Dev. 1977, 16, 157. (37) Flow of Fluids through Valves, Fittings, and Pipe; Technical Paper No. 410; Crane Company: Stamford, CT, 1988. (38) Kannuluik, W. G.; Carman, E. H. The Temperature Dependence of the Thermal Conductivity of Air. Aust. J. Sci. Res. 1953, 4, 305.

3294


Multiobjective Optimization of a Fixed Bed Maleic Anhydride Reactor

Recommend Documents