5196
Ind. Eng. Chem. Res. 1997, 36, 5196-5206
Optimization of the Reaction Conditions for Complex Kinetics in a Semibatch Reactor Juha Lehtonen,† Tapio Salmi,*,† Antti Vuori,‡ and Heikki Haario§ Laboratory of Industrial Chemistry, A° bo Akademi, FIN-20500 Turku, Finland, Kemira Agro Oy, Research Centre, FIN-02271 Espoo, Finland, and Department of Mathematics, University of Helsinki, FIN-00100 Helsinki, Finland
Optimization of the reaction conditions in liquid phase semibatch reactors where several complex reactions simultaneously take place was studied. The optimization procedure was demonstrated with the production of substituted alkyl phenols from aromatic sulfonic acid via alkali fusion. The reagent, the substituted aromatic sulfonic acid, participates in the main reaction and in several side reactions, such as the elimination of the alkyl, the amino, and the sulfonic acid groups as well as polymerization reactions. The kinetics was determined for the reaction system, on the basis of the laboratory scale batch reactor experiments. The kinetic model was combined with a large scale reactor model, where not only the feed of the reactant but also the precipitation and volatilization of some reacting compounds were considered. It turned out that an extrapolation from the 0.1 dm3 laboratory reactor to the 10 m3 large scale was very successful. The reaction conditions of industrial scale semibatch reactor were optimized using the desirability function approach, which combines the desired conversion of the reagent, the desired selectivity of the product, and the reaction time. The qualities of different desirability functions were compared. The great benefit of the desirability approach is that it provides smooth objective functions, avoiding a cumbersome constrained optimization. Introduction Fine and speciality chemicals are often produced in semibatch reactors. The reaction system in this kind of small scale production is often very complex, including several consecutive-competitive reactions. Optimization of these systems is a difficult task; there exists no general rules for the optimization of reaction conditions in complex reaction networks. In this study, unconstrained optimization is applied on a complex reaction network, the network of the alkali fusion of a polysubstituted aromatic ring. Unconstrained optimization has many advantages compared to the traditional constrained optimization, e.g., serious numerical problems can be avoided in this way. The optimization will in the current work be carried out with desirability functions, the definitions of which are based on desired conversions and selectivities. Aromatic hydroxy compounds are important intermediates in the dye-stuff industry (Booth, 1988). A typical synthesis route for the production of hydroxy compounds proceeds via alkali fusion, where the reagent, the sulfonated aromatic compound, undergoes a substitution reaction with alkali metal (M) hydroxide. The reaction takes place in the liquid phase, in the presence of an excess of alkali metal hydroxide. For instance, the reaction of a substituted aminosulfonic acid with alkali is described as follows:
The alkali metals present are sodium (Na) and potassium (K). Potassium is used as an additive to †
A° bo Akademi. Kemira, Research Centre. § University of Helsinki. ‡
S0888-5885(97)00244-3 CCC: $14.00
suppress the melting point of the hydroxide: a mixture containing 39 mol % NaOH and 61 mol % KOH has an eutectic point at 185 °C, whereas the melting point of pure NaOH is 360 °C. The reaction also produces alkali metal sulfite, which precipitates as a solid phase. Alkali fusion is principally a very old industrial process, being applied to both mono- and diaromatic sulfonates as reviewed by Booth (1988). Typical monoaromatic intermediates produced through alkali fusion are maminophenol and N,N-diethyl-m-aminophenol. Alkali fusion is industrially carried out at temperatures of 200-400 °C in a mixture of NaOH and KOH (Wedemeyer, 1976). Because of safety aspects, i.e., the exothermic reactions, the large scale processes are operated in semibatch mode: the aromatic reagent and water are fed into the alkali melt. The industrial process is complicated by the presence of side reactions, which partially destroy the aromatic reagent and partially convert the desired product to less valuable byproducts. Typical side reactions are the eliminations of the sulfonic, the alkyl (R and R′), and the amino groups as well as the oligomerization of the aromatic rings. The crucial parameters to be optimized in a semibatch reactor are the reaction time, the feed time, and the reaction temperature. The endeavor of this work is to develop the optimization procedure of a complex reaction network in semibatch conditions by using the alkali fusion as a demonstration example. The reaction system was analyzed using a reaction scheme and a rigorous kinetic model based on experimental data from laboratory batch reactor experiments. A complete kinetic model for the alkali fusion has been presented by Lehtonen et al. (1995). Mechanisms and Kinetics of Alkali Fusion The principal reactions in the alkali fusion process are as follows: (1) substitution of the sulfonic acid group with a hydroxyl group, (2) elimination of the sulfonic © 1997 American Chemical Society
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5197
acid group, (3) dealkylation of R′, (4) dealkylation of R, usually -CH3 (5) elimination of the amino group, and (6) polymerization reactions. The main reaction (1) can be written as
(1)
The elimination and the dealkylation reactions (25) are illustrated below:
R′3 ) k′3cA
(9)
The dealkylation of R′ (eq3) is represented by reactions 9-14 in Figure 1. The dealkylation of the alkyl group R (eq 4) has been proved to proceed through the formation of a carboxylate ion as an intermediate: for the case that R was a CH3 group, traces of substituted aromatics with COOH groups were detected experimentally. The amounts of the intermediates were, however, almost negligible in all experiments, which suggests that they react very rapidly to the final products. Thus it was sufficient to describe the kinetics of reactions 15-21 in Figure 1 with a simple second-order rate law:
R′4 ) k′4cAcOH
(10)
The elimination of the amino groupsreaction category 5swas observed only for the phenolic compounds as illustrated in Figure 1. Thus this elimination contributes in four reactions, 22-25, in the scheme. These reactions were assumed to be monomolecular thermal decompositions having the kinetics
R′5 ) k′5cA
Several polymerization reactions (formation of tar) can take place in the reaction milieu; the most important polymerization process in the actual reaction mixture is, however, a bimolecular reaction between two aromatic rings:
(11)
The formation of tar is a polymerization process principally taking place with all aromatic rings leading to dimers, trimers, etc. For the present case it turned out, however, that it is sufficient to consider the dimer formation and ignore the formation of trimers, tetramers, etc. The dimerization process was supposed to be bimolecular, and all the dimerization rate constants were presumed to be equal. The rate law becomes
R′6 ) k′6
∑ ∑
ckc1
(12)
k)A...N 1)k...N
The reaction routes of the aromatic rings are illustrated in Figure 1. For the mathematical modeling of the large scale reactor it is necessary to know the kinetics of the main and the side reactions (1-6). The kinetic experiments carried out isothermally at 220270 °C in closed batch reactors revealed that the scheme depicted in Figure 1 is suitable for the description of the reactions. All reactions can in the actual conditions be considered to be irreversible. The main reaction (1) is a bimolecular nucleophilic substitution; thus a secondorder rate law was applied
R′1 ) k′1cAcOH
(7)
where cA denotes the concentration of any aromatic ring undergoing the substitution (reactions 1-4 in Figure 1). Also the elimination of the sulfonic acid group, reaction (2), proceeds via a bimolecular process requiring the presence of hydroxide ions. Consequently, secondorder kinetics was applied on these reactions:
R′2 ) k′2cAcOH
A detailed treatment of the reaction scheme depicted in Figure 1 would require the use of 25 different rate constants for reactions 1-25. This is far too much, regarding the amount of experimental data available. Thus reasonable simplification is to apply the principle of independent reactivity of a functional group: the reactivities of the -SO3H, -NHR′, -R, and -N groups were assumed to be independent of the remaining aromatic ring. This is principally an oversimplification, but it is reasonable for the present case. It should be recalled that analogous simplifications have been used for years in polymerization kinetics. The primary rate constants for reactions 1-25 (k1-k25) were thus grouped according to the principal reaction categories 1-6:
(8)
The elimination reactions are represented by steps 5-8 in Figure 1. The dealkylation of R′ was regarded as a thermal decomposition, which suggests the use of a first-order rate law:
k′1 ) k1 ) k2 ) k3 ) k4
(13)
k′2 ) k5 ) k6 ) k7 ) k8
(14)
k′3 ) k9 ) k10 ) k11 ) k12 ) k13 ) k14
(15)
k′4 ) k15 ) k16 ) k18 ) k19 ) k20 ) k21
(16)
k′5 ) k22 ) k23 ) k24 ) k25
(17)
k′6 ) k27
(18)
Experimental observations revealed that the polymerization process played the most important role in the beginning of alkali fusion. This implies that polymerization proceeds mainly between the aromatic com-
5198 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997
Figure 1. The reaction scheme for the alkali fusion (M ) Na or K, R and R′ are alkyl groups).
pounds A and B being present in the melt phase at the initial stage of the reaction. Thus the reaction kinetics (12) consists of three contributions (A-A, A-B, B-B) only. After these simplifications the rates of reactions 1-28 can be written as follows: R1 ) k1cAcO R2 ) k2cBcO R3 ) k3cCcO R4 ) k4cDcO R5 ) k5cAcO R6 ) k6cBcO R7 ) k7cCcO R8 ) k8cDcO R9 ) k9cA R10 ) k10cE R11 ) k11cI R12 ) k12cC R13 ) k13cG
(19) (20) (21) (22) (23) (24) (25) (26) (27) (28) (29) (30) (31)
R14 ) k14cK R15 ) k15cAcO R16 ) k16cEcO R17 ) k17cIcO R18 ) k18cBcO R19 ) k19cFcO R20 ) k20cJcO R21 ) k21cMcO R22 ) k22cE R23 ) k23cF R24 ) k24cG R25 ) k25cH
(32) (33) (34) (35) (36) (37) (38) (39) (40) (41) (42) (43)
The rate equations for the dimerization reactions (6) are:
R26 ) k26cA2
(44)
R27 ) k26cB2
(45)
R28 ) k26cAcB
(46)
The generation rates of the compounds (ri) were obtained from the reaction rates (Rj) using the stoichiometry s
ri )
νijRj ∑ j)1
(47)
where νij is the stoichiometric coefficient of component i in reaction j. The stoichiometric coefficients originated from the reaction scheme (Figure 1) are given in Table 1. The numerical values of the rate constants including the temperature dependences were determined with
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5199 Table 1. The Stoichiometric Matrix
parameter estimation are shown in Table 2. The results were presented more completely by Lehtonen et al. (1995). Experiments in Industrial Scale
nonlinear regression analysis (MODEST software; Haario, 1994) from the batch reactor data. Determination of the Kinetic Parameters in Alkali Fusion The kinetic experiments were carried out in a laboratory scale pressurized autoclave of Hastelloy B (0.100 dm3). The experiments were performed at 220, 230, 240, 250, 260, and 270 °C. The main analytical method for obtaining the kinetic data was high-pressure liquid chromatography (HPLC). The experimental setup and the analytical methods are described in detail in Lehtonen et al. (1995). An example of the experimental data from the batch reactor and the fit of the kinetic model are given in Figure 2, where the concentrations of A, B, E, F, and S at 240 and 260 °C are depicted. As can be seen from the figure , the kinetic model is able to reproduce the main trends of the experimental data; for example, the concentration maximum of the main product, compound E. Figure 2 shows convincingly that even the simplification concerning the independent reactivity of a functional group is very justified in the present case (eqs 13-18). Without this kind of simplification the system would be terribly overparametrized and a meaningful parameter estimation would be impossible. In order to further investigate the reliability of the kinetic model, the temperature dependences of the kinetic parameters were checked by preparing the Arrhenius’ plots. The Arrhenius’ plots of parameters k′1, k′2, k′3, k′5, and k′6 are depicted in Figure 3. As can be seen from the plot, the rate parameters follow the law of Arrhenius reasonably well. The agreement for parameters k′1 (the main reaction), k′2 (elimination of -SO3H), and k′6 (tar formation) were excellent, whereas the temperature dependences of parameters k′3 (dealkylation of R′) and k′5 (elimination of NH2) were slightly obscured by experimental scattering. The role of reaction 4 (dealkylation of R) turned out to be very minor, so this reaction was discarded. The final values of the rate constants were obtained by merging all data sets together and applying nonlinear regression analysis directly on the primary data. The final results of the
Three experiments were carried out in an industrial scale stirred-tank reactor with a volume of 10 m3. The reaction vessel was made of nickel. The reactor temperature was controlled with an oil jacket. The reactor was operated in semibatch mode. Alkali hydroxide (MOH) and water were charged into the reactor prior to the experiment, the heating was started, and after attaining the desired temperature, the reagent (A) was fed into the reactor as an aqueous solution. Some water was continuously removed from the liquid phase and the alkalisulfite was precipitated, but it remained in the reaction vessel. Samples were withdrawn from the reactor after the feeding period. The temperature measurement was considered to be quite reliable after the feeding period, but during the period the cold feed induced hightemperature gradients in the reactor and it was very difficult to obtain reliable information on the reactor temperature. The correction for the feed temperature as well as for the temperature after the feed period was estimated with the kinetic model. It was found out for all three experiments that the estimated temperature was 22 °C lower during the feed and 8 °C lower after the feed than the measured temperature. The experimental conditions for the industrial scale experiments I, II, and III are reported in Table 3. The simulated and the experimental concentrations are presented in Figure 4, which shows that the model based on a 0.1 dm3 scale agrees with the experimental concentrations of a 10 m3 scale surprisingly well. Mass Balances for a Semibatch Reactor The general mass balance for an arbitrary compound (i) in the liquid phase of a semibatch reactor can be written as
rim )
dni + dt
∑n˘ ′i
(48)
where ri is the generation rate of i per unit mass of reaction liquid, m is the mass of the liquid phase, ∑n˘ ′i is the sum of the flows into the liquid phase and out from the liquid phase (n˘ ′i < 0 for inlet flows and n˘ ′i > 0 for outlet flows). In the present case the equation becomes
n˘ i,in + rim ) n˘ i,out + n˘ ′i,out +
dni dt
(49)
where n˘ i,in denotes the feed rate of a compound i, n˘ i,out is the flow out to the gas phase due to volatilization, and n˘ ′i,out is the amount of substance per time precipitated in the sulfite phase. The mass-based concentrations are introduced:
ni ) cim
(50)
n˘ ′i,out ) c′iΦ
(51)
and
5200 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997
Figure 2. The concentration curves of components A, B, E, F, and S at 240 °C (a) and 260 °C (b).
where Φ is the mass precipitated per unit of time. Consequently, the accumulation term can be expressed as
β)
mR m0
(55)
dni dci dm ) m + ci dt dt dt
γ)
mD m0
(56)
δ)
mQ m0
(57)
(52)
For the reagent (mR), the reactor (m), the precipitated (mQ), and the distillate (mD) masses, the following total balance is valid:
mR + m + mQ + mD ) mR,0 + m0
(53)
where the terms in the right-hand side describe the initial state of the system. To simplify the subsequent treatment, some dimensionless quantities are introduced:
R)
mR,0 m0
(54)
After introducing definitions (54)-(57) in eq 53, a relation between the reactor mass and the initial mass is obtained:
m )1+R-β-γ-δ m0
(58)
Correspondingly, the mass flows can be expressed with β, γ, and δ:
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5201
dmR dβ ) -m0 dt dt
Table 2. Rate Constants
(59)
m ˘ in ) m ˘ out )
dmD dγ ) -m0 dt dt
(60)
dδ dt
(61)
Φ ) m0
rate constant
A
Ea
k′1 k′2 k′3 k′5 k′6 k17
2.47 × 1014 kg/mol‚h 3.96 × 107 kg/mol‚h 9.71 × 1016 1/h 1.17 × 101 1/h 6.63 × 1020 kg/mol‚h 3.16 × 109 kg/mol‚h
166 kJ/mol 103 kJ/mol 195 kJ/mol 40.1 kJ/mol 225 kJ/mol 102 kJ/mol
After introducing the definitions (54)-(61), the original mass balance equation (49) becomes
dβ dγ + (1 + R - β - γ - δ)ri ) ci,out + dt dt dci dδ dβ dγ dδ c′i - ci + + + (1 + R - β - γ - δ) dt dt dt dt dt (62)
-ci,in
(
)
For a further development of eq 62 it is practical to consider the different categories of the compounds separately. For the reagent molecules (A) fed into the reactor we have c′A ) 0 and cA,out ) 0. The solution of dcA/dt from (62) gives
dcA (cA - cA,in) dβ ) rA + + dt 1 + R - β - γ - δ dt cA dγ dδ + (63) 1 + R - β - γ - δ dt dt
(
)
For the other compounds, B-O and R-S, c′i ) 0, ci,out ) 0, and ci,in ) 0, and the mass balance equation (62) becomes
ci dci dβ dγ dδ ) ri + + + dt 1 + R - β - γ - δ dt dt dt
(
)
(64)
For water (P) c′P ) 0, and the mass balance obtains the form
dcP (cP + cP,in) dβ ) rP + + dt 1 + R - β - γ - δ dt (cP - cP,out) dγ dδ + cP (65) 1 + R - β - γ - δ dt dt
Figure 3. The Arrhenius’ plot of the rate constants k′1, k′2, k′3, k′5, and k′6. Table 3. The Conditions of the Industrial Scale Experiments (dmR/dt)/ (kg/h)
T/°C
0-19.25 19.25-23.75 23.75-34.75 34.75-36.75 36.75-44.25 44.25-45.25 45.25-49.25 49.25-50.25 50.25-51.25
218 218 f 250 250 250 f 255 255 255 f 240 240 240 f 250 250
Experiment I 484.46 51.25-52.25 0 52.25-68.25 0 68.25-70.75 0 70.75-73.75 0 73.75-74.25 0 74.25-74.75 0 74.74-76.25 0 76.25-86.00 0
228 228 f 250
Experiment II 777.53 24-28 0 28-36.25
228 228 f 250 250 f 260 260 f 264
Experiment III 817.39 31-35 0 35-37 0 37-39.5 0
0-9.5 9.5-24 0-9 9-24 24-27 27-31
T/°C
(dmR/dt)/ (kg/h)
250 f 255 255 255 f 260 260 260 f 265 265 f 260 260 260 f 250
0 0 0 0 0 0 0 0
250 f 265 265 f 253
0 0
264 f 266 266 266 f 265
0 0 0
t/h
t/h
For sulfite (Q) cQ,in ) 0 and cQ,out ) 0, which gives
dcQ cQ dβ dγ dδ ) rQ + + + dt dt dt (1 + R - β - γ - δ) dt
(
)
(66)
On the other hand, the sulfite balance is used to obtain an expression for dδ/dt. Multiplication of the original balance equation (49) (n˘ Q,in ) n˘ Q,out ) n˘ ′Q,out ) 0) with the molar mass of sulfite (MQ) gives (nQMQ ) mQ)
dmQ ) rQmMQ dt
(67)
For water the pseudo-steady-state approximation based on the requirement of thermal balance of the reactor is applied: the water fed into the reactor and the water generated in the chemical reaction are set equal to the amount of water leaving the reactor; an accumulation would cause a dramatic temperature decrease:
m ˘ incP,in + rPm ) n˘ P,out
(69)
The mass flow of water (m ˘ P) is solved from (69) after multiplication with the molar mass:
which after division with the original mass (m0) becomes
dδ ) rQ(1 + R - β - γ - δ)MQ dt
˘ incP,in + rPm)MP m ˘ P ) (m
(70)
(68)
i.e., the precipitation of the sulfite is assumed to proceed with an infinite speed: all the sulfite, which is generated, is immediately converted into the solid phase.
On the other hand, we have the following definitions:
m ˘ P)
dmD dt
(71)
5202 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997
The mathematical model of the large scale reactor consists of the mass balance equations of the compounds, (63)-(66), and the differential equations for δ and γ, eqs (68) and (73). The following initial conditions were applied at t ) 0:
ci(0) ) c0,i, i ) A...S
(74)
γ(0) ) 0
(75)
δ(0) ) 0
(76)
Objective of the Optimization In the optimization of the reaction conditions several criteria can be used. Basically, in industrial applications the criterion is the profitability of the whole process, including usually several reaction and separation steps. When the reactions are carried out in a semibatch reactor, the crucial optimization parameter affecting the profitability is the total time spent on yielding the desired product quality (productivity). The total time is mainly affected by the reaction time, the reaction temperature, and the preparation time between the batches. In a multicriterium situation, the overall optimum usually is a compromise between contradicting objectives. There might indeed be no objectively “best” solution, but the investigator must weigh the different factors in a more or less subjective way. Here, the desirability function technique was employed. The idea is to specify separately the “good” and “bad” values for each target by, e.g., graphical means, and then pool them together into an overall desirability function. As discussed in the next section, various functions (eqs 81 and 82) may be used for the desirability profiles which convert the response (S) into the desirability value (D), 0 < D(S) < 1. Here 0 indicates a completely unsatisfactory result, while 1 indicates that the required level of a response has been reached. The overall desirability is then defined as the product (or geometric mean) of the individual target profiles. A great variety of targets and constraints may be formulated as desirabilities, and a complex optimization with general constraints may be performed as a straightforward optimization with simple bounds. In the present case, the objective function to be maximized is given by the following variant of the idea:
f)
Figure 4. The industrial scale experiments I (a), II (b), and III (c) in the conditions given in Table 3.
m ˘ in )
dmR dt
(72)
Division of eq 70 with m0, recalling the definitions (54)(57) and the relations (71)-(72) gives finally a relationship between γ and β:
dβ dγ ) - cP,in + rP(1 + R - β - γ - δ) MP dt dt
(
)
(73)
The amount of the feed solution (β) was measured experimentally. This pseudo-steady-state approximation for the volatiles (water) can only be applied for this particular case. A general model for a distillation process connected to a semibatch reactor is presented by Lehtonen et al. (1997).
DXADSE tr + tD
(77)
where DXA denotes the desirability of the conversion of the reactant A, DSE denotes the desirability of the selectivity of the desired product E, tr is the reaction time, and tD is the preparation time. The conversion and the selectivity are defined accordingly:
n0A - nA cA )1(1 + R - β - γ - δ) n0A c0A
(78)
nE cE(1 + R - β - γ - δ) ) n0A - nA c0A - cE(1 + R - β - γ - δ)
(79)
XA )
SE )
The product of the conversion and the selectivity gives the yield of the desired compound,
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5203
YE )
nE ) XASE n0A
In the present case, the main interest should be focused on the upper left corner of the reaction scheme (Figure 1), since the dominant reactions take place between compounds A, B, E, F, and S (tar). These reactions form the following simplified scheme: k′1
k′3
(I) A 98 E 98 F k′3
D)
(80)
k′1
(II) A 98 B 98 F k′6
(III) A 98 S In this scheme, there exist two competitive-consecutive reaction pathways to F, where the reaction intermediate of the reaction (I) E, is the desired product. As a competitive process reaction (III) forms tar. There might exist an optimal temperature, which maximizes the yield of the desired product (E). According to Table 2, the rate constant k′6 has the highest activation energy. Also the activation energy of k′3 is higher than the activation energy of the desired main reaction constant (k′1). Although the difference between the activation energies is not very high, it is assumed that a lower temperature will result in a better selectivity, due to a lower tar formation rate and lower decomposition rates of E. On the other hand, higher temperature will give a higher reactant conversion in a shorter reaction time. Optimal temperature schedules for batch reactors based on the Pontryagin’s principle were given, for example, by Rippin (1978). These useful recommendations are valid for simple reaction systems and batch reactors. In this kind of complex reaction system carried out in a semibatch reactor, it is very hard to try predict optimal temperature profiles with the methods of exact mathematical analysis. Selection of Desirability Functions The optimization problem is typically a constrained one, i.e., the conversion and selectivity should remain within some limiting values. In practice, these quantities are continuous functions of the reaction conditions, e.g., the temperature. Therefore, the limits of these quantities are de facto smooth, and a rigorous choice of the constraints inevitably becomes a subjective procedure. Furthermore, constrained optimization often suffers from serious numerical problems. Consequently, a more flexible alternative to the traditional constrained optimization becomes attractive. One possible way to avoid the constrained optimization is to introduce the concept of desirability functions, where continuous desirabilities are defined for these quantities. For example, trigonometric and exponential functions can be used as desirability functions. Of course, the conversion and the selectivity can also be optimized directly, but very unwanted results can appear, e.g., a reaction temperature or a reaction time which cannot be realized in practice. The choice of desirability functions can easily become a subjective procedure. In this paper some criteria for the choice of the desirability functions are presented. The following exponential and trigonometric functions are proposed for desirability functions:
D)
1 1 + e-(S-S0)/R
1 arctan(β(S - S0)) + 1/2 πR
(81) (82)
where S0 is the value (e.g., conversion, yield, or selectivity) which gives the desirability 0.5. Parameters R and β define the slope of the curve at S0 (Figure 5). The slope of function (81) is given by 1/4β and the slope of function (82) is Rβ/π. In practice, just two parameters should be defined: S0 and the slope. When functions like (82) are used, the values of the functions at S ) 0.0 and at S ) 1.0 are fixed with parameter R. In order to fulfill the criteria D e 0.01, when S ) 0.0 (R1) and D g 0.99 when S ) 1.0 (R2) the following expressions are introduced:
arctan(-β1S0) -0.49π
(83)
arctan(β2 - β2S0) -0.49π
(84)
R1 e and
R2 e
and R1 ) R2 only if S0 ) 0.5. If the value S0 * 0.5 is given, the function should consist of two parts: the first part ranging from 0.0 to S0 and the second part from S0 to 1.0. The derivatives (slopes) of both parts at S0 should coincide. The behavior of function (81) at S ) 0.0 and S ) 1.0 is somewhat better than the behavior of function (82), since it approaches steeper 0 and 1. Anyway, in some cases the values of the function at these points differ decisively from 0 and 1. In these cases it is recommendable to use functions as (82). Optimization of the Alkali Fusion In our case study, the exponential as well as the trigonometric functions with the following parameters were introduced as desirabilities for the conversion and selectivity: DXA:
exponential, S0 ) 0.925, β ) 0.005 trigonometric, S0 ) 0.925, R ) 1.016, β ) 154.626 when XA e 0.925 S0 ) 0.925, R ) 0.967, β ) 162.461 when XA > 0.925 DSE:
exponential, S0 ) 0.5, β ) 0.05 trigonometric, S0 ) 0.5, R ) 0.943, β ) 16.660 The forms of these desirability functions are illustrated in Figure 5. The numerical parameters in the equations given above were adjusted in such a way that a desirability of 0.5 was achieved with the conversion XA ) 0.925 and with the selectivity SE ) 0.5. The objective function for the optimization was defined by (77).
5204 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 Table 4. The Optimization Results the isothermal optimization the optimized temperature profile isothermal 260 °C isothermal 295 °C trigonometric desirability functions optimization without desirability functions
tr/h
tfeed/h
Tr/ °C
Tfeed/°C
XA
SE
f
DXA
DSE
7.91 7.33 7.33 7.33 8.27 2.48
6.89 6.33 6.33 6.33 7.19 2.22
283.3 profile 260 295 282.6 301.4
282.7 profile 260 295 282.0 299.5
0.943 0.944 0.654 0.983 0.943 0.940
0.482 0.485 0.456 0.411 0.483 0.430
0.234 0.253 5.16 × 10-25 0.088 0.214
0.972 0.977 2.90 × 10-24 0.999 0.908
0.409 0.424 0.292 0.144 0.408
Figure 6. The value of the objective function as a function of the reaction temperature (a) and the reaction time (b). Figure 5. The desirability functions DXA (a) and DSE (b) by using exponential and trigonometric functions.
The definition of the objective function implies that f increases monotonically with XA and SE, and the product of XA and SE gives the yield according to eq 80. Furthermore, the objective function decreases with a longer reaction time. In order to avoid the weighing of the reaction time, the sum tr + tD was normalized to vary between zero and unity. The objective function was maximized by using the software MODEST (Haario, 1994), where a modified Simplex algorithm (Nelder and Mead, 1965) was used for the optimization and the differential equations were solved with the LSODE subroutine (Hindmarsh, 1983). In the first optimization, a constant temperature for the feed period and for the period after the feed as well as the feed and the total reaction time were optimized. The results are presented in Table 4. The sensitivities of the objective function as a function of the temperature and the reaction time are presented in Figure 6. As can be seen from the figure, both quantities are very sensitive to the changes in the temperature and the reaction time. On the basis of the result of the first optimization, the feed time was fixed to tr - 1 h and the second optimization was performed. In this case the
Figure 7. The optimized temperature profile.
reaction time was divided into seven steps. The temperatures of these steps were optimized as well as the reaction time. The temperature between the steps was interpolated linearly. The optimized temperature profile is presented in Figure 7. According to this profile, a quite high temperature is preferable in the beginning, after which the temperature should be decreased until the feed is commenced. After the feed period, the
Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 5205
Conclusions
Figure 8. The simulated concentrations of components A (a) and E (b) at 260 °C, 295 °C, and the optimized temperature profile.
temperature should be raised again. The optimized concentrations as well as simulated isothermal concentrations of the reactant A and the product E at 260 and 295 °C are presented in Figure 8. The temperature decrease during the feed suppresses the major side (tar formation) reactions, because these reactions have a higher activation energy than that of the main reaction (Table 2). On the other hand, the temperature can be increased after the feed period (Figure 7), since the roles of the tar formation reactions become less important because of the consumption of the reagent (A). In the first optimizations, exponential desirability functions were used. In order to test the trigonometric desirability functions and the optimization without desirability functions, two further optimizations were carried out. These optimizations were performed analogously with the first optimization with constant feed and reaction temperatures. The results are presented in Table 4. The difference between using the exponential and the trigonometric desirability functions is not remarkable. By using the trigonometric function, a slightly longer reaction time was proposed. On the other hand, the optimization without desirability functions gave a totally different result, with a lower selectivity and conversion. Also the proposed reaction temperature and the reaction time are unrealistic from the thermodynamic point of view. In the discussion of the temperature optimization it is important to compare the optimal temperature profile with the simplest alternative, an isothermal batch. According to Table 4 the optimal temperature profile guarantees the highest selectivity of E. Of course a higher temperature gives a higher conversion and a shorter reaction time, but due to the proper selection of the desirability function, also the need for a high selectivity can be taken into account.
The reaction conditions of a complex organic reaction system in an industrial scale semibatch reactor were optimized. The production of substituted alkylphenols using alkali fusion was chosen for a demonstration system (Figure 1). A kinetic model was developed for the reactions taking place in the alkali fusion. The analysis of the laboratory scale batch reactor data revealed that the principle of constant reactivity of a functional group is applicable. Using the kinetic model and the model for the industrial scale semibatch reactor, three experiments were simulated and a very good argeement with the experimental results (Table 3 and Figure 4) was obtained. The optimization procedure was carried out using the unconstrained optimization and desirability functions. A procedure for the choice of the suitable desirability functions was introduced and different desirability functions were compared in the optimization (Figure 5). The optimization was proceeded stepwise, obtaining first an isothermal optimum and optimal reaction and feed time. After that an optimal temperature profile was calculated for the system (Figure 6). These results were compared with some isothermal runs as well as the optimization without desirability functions (Table 4 and Figure 7). It can also be concluded that the unconstrained optimization without desirability functions usually gives results, which are unrealistic in practice. This paper presents a proper procedure to optimize a semibatch reactor. The startpoint in this study is the kinetic model, and the optimization is mainly based on the model developed in the laboratory scale. Although this may be the best way to solve the scale-upoptimization problem, there may other points, e.g., the reaction thermodynamics, which should be considered before the full scale production is started. Acknowledgment The skillful help of Ms. Paula Nousiainen in the experimental work is gratefully acknowledged. Notation A ) frequency factor c ) concentration (mol/kg) D ) desirability Ea ) activation energy f ) objective function, eq 77 k, k′ ) rate constant m ) mass m ˘ ) mass flow M ) molar mass n ) amount of substance n˘ , n˘ ′ ) flow of substance r ) generation rate of a compound Q ) objective function R, R′ ) reaction rate S ) selectivity T ) temperature t ) time (tr ) reaction time, tD ) preparation time) X ) conversion Y ) yield Greek Letters R, β, γ, δ ) dimensionsless masses R, β ) coefficients in the desirability functions ν ) stoichiometric coefficient ν ) stoichiometric matrix Φ ) mass precipitated per time
5206 Ind. Eng. Chem. Res., Vol. 36, No. 12, 1997 Subscripts and Superscripts 0 ) initial value D ) distillate i ) component index j ) reaction index in ) incoming OH ) hydroxide R ) feed solution out ) outgoing Abbreviations A, B, ..., S ) reaction components M ) alkali metal ion R ) alkyl group R′ ) alkyl group
Literature Cited Booth, G. The Manufacture of Organic Colorants and Intermediates; Society of Dyers and Colourists: Bradford, 1988; pp 3167. Haario, H. MODEST-User’s Guide; Profmath Oy: Helsinki, 1994. Hindmarsh, A. C. ODEPACK-A Systematized Collection of ODESolvers. Scientific Computing; Stepleman, R., et al., Eds.; IMACS/North-Holland Publishing: Amsterdam, The Netherlands, 1983; p 55-64.
Lehtonen, J.; Salmi, T.; Vuori, A.; Haario, H.; Nousiainen, P. Modeling of the Kinetics of Alkali Fusion. Ind. Eng. Chem. Res. 1995, 34, 3678-3687. Lehtonen, J.; Salmi, T.; Harju, T.; Immonen, K.; Paatero, E.; Nyholm, P. Dynamic Modelling of Simultaneous Reaction and Distillation in a Semibatch Reactor System. Chem. Eng. Sci. 1997, in press. Nelder, J. A.; Mead, R. A Simplex method for function minimization. Comput. J. 1965, 7, 308-13. Rippin, D. W. T. Improvements Obtained in the Performance of Batch Reactors by the Adjustment of Operating Conditions during the Course of the Reaction. SIA/FVC Tagung, Basel, Switzerland, 1978 (SEG/V/13/78, TCL-ETH, Zu¨rich). Wedemeyer, K.-F. Die Sulfonsa¨ure-Gruppe. Methoden der Organishen Chemie, 4th ed.; Mu¨ller, E., Ed.; Georg Thieme Verlag: Stuttgart, 1976; pp 202-229.
Received for review March 19, 1997 Revised manuscript received June 27, 1997 Accepted July 12, 1997X IE9702441
X Abstract published in Advance ACS Abstracts, October 1, 1997.
