A systematic generation of reactor designs

Computers and Chemical Engineering 29 (2005) 1101–1112

A systematic generation of reactor designs II. Non-isothermal conditions Magne Hillestad∗ Cybernetica AS, Leirfossvegen 27, N-7038 Trondheim, Norway Received 22 January 2004; received in revised form 19 November 2004; accepted 30 November 2004

Abstract A new method is proposed for systematic generation of conceptual design of reactor networks. Here, non-isothermal conditions are considered. Optimal distribution of heat transfer area, in addition to optimal mixing and distribution of feeds are targeted. Given temperature and composition of feed streams, the objective is to find the sequence and size of ideal reactors, the distribution of extra feed streams, the heat transfer distribution along the reactor path and the total reaction time that maximizes the space–time yield of the key product component. The method can be applied to solve problems of any number of components and reactions. The method is applied to industrially important processes such as the methanol synthesis and the steam methane reforming process. Interesting results are obtained with respect to the required heat exchange area. © 2004 Elsevier Ltd. All rights reserved. Keywords: Optimal reactor design; Reactor network synthesis; The maximum principle; Optimal mixing; Heat transfer; Distributed feed; Methanol synthesis; Steam methane reforming

1. Introduction Given a chemical system of reactions, a description of their kinetics and given feed composition and temperature, how should the reactants be contacted or mixed so as to optimize a production objective. This is still one of the great problems of chemical reaction engineering. The reactor design or synthesis involves determining the type, size and interconnections of reactor units and distribution of feeds. In addition, the heat transfer is an important issue for most industrial chemical reactors. The design activity can greatly impact the process energy use, separation requirements and overall economics. As waste minimization and sustainable process technologies are increasingly becoming more important, more focus on the reactor design problem is required. The real challenge of this problem arises when complex reaction schemes are encountered, such as parallel and consecutive reactions. Adding to the complexity, is the fact that the reactions taking place affect the reaction tempera∗

Tel.: +47 73 82 28 76; fax: +47 73 82 28 71. E-mail address: [email protected].

0098-1354/$ – see front matter © 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2004.11.009

ture, which again affect the reaction rates. A reactor design method should be able to handle complex reactions and rate– temperature interactions. In addition to the evolutionary modification normally applied in the process industry, the efforts in academic research can be classified into two general categories; superstructure based optimization and systematic generation methods. A superstructure is a general process network that incorporates a predefined number of units and their interconnections. It gives rise to a formulation that is normally solved by mixed integer and nonlinear programming. Early attempts include Jackson (1968) who employed optimal control techniques to optimize a reactor network superstructure composed of plug flow reactors connected by side streams. Belonging to the systematic generation methods is the attainable region method (AR), on which a large number of articles have been published. A recent article by Feinberg (2002) gives an outlook of the attainable region method applied on pure reactor synthesis and reactor–separator synthesis. Another AR article dealing with non-isothermal condition is given by Hopley, Glasser, and Hildebrandt (1996). Balakrishna and Biegler (1992) formulated a targeting model solved by nonlinear programming.

1102

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

Nomenclature a A AR A∗ c CSTR ek E f h H J L m nc nx N NSt PFR r r R ˜ R RT Rk STY T TL Tw T uF uH uH uM uT U

specific heat transfer area, a = AR /VR [m2 /m3 ] heat transfer area increasing with space–time [m2 ] total heat transfer area [m2 ] relative heat transfer area, A∗ = (UA)/(UA)R [–] concentration vector of dimension (nc × 1) [mol m−3 ] continuously stirred tank reactor the kth unit vector, e.g. e2 = [0, 1, 0, . . . , 0] a matrix of dimension (nx × nx ), E = diag(0 · · · 0, 1) the model formulation comprising both PFR and CSTR defined by (Eq. 17), f = dx/dξ vector of reduced heat of reaction, h = (−H)/( cp ) [K m3 mol−1 ] vector of heat of reaction [J mol−1 ] the Jacobian matrix, J = ∂R/∂x the partial derivative, L = ∂f/∂x number of reactions, m = dim(r) number of components, nc = dim(c) dimension of x, nx = dim(x), nx = nc + 1 stoichiometric matrix augmented with heat of reactions, dimension (nx × m) Stantons number for heat transfer, NSt = τR /τH [–] plug flow reactor number of reactors along the path vector of reaction rates [mol m−3 s−1 ] vector of component formation rates augmented with rate of heat formation, R = Nr ˜ x) while R is the average reaction rate, R(θ, includes deactivation dynamics rate of heat formation, RT = h r [K s−1 ] formation rate of the key component k [mol m−3 s−1 ] space–time yield of the key component reactor temperature [K] time between catalyst reload [s] temperature of cooling/heating medium [K] temperature difference T = T − Tw [K] the design function defining the distribution of feeds [–] the design function for heat transfer [s−1 ] vector containing the parameterized design function uH the design function defining the mixing regime [–] the design temperature of the cooling/heating medium [K] overall heat transfer coefficient [Wm−2 K−1 ]

U VR x y

diagonal matrix of dimension (nx × nx ), U = diag(0 · · · 0, NSt ) reaction volume [m3 ] vector of concentrations and temperature, dimension (nx × 1), x = [y , T ] vector of mole fractions, dimension (nc × 1)

Greek symbols θ ξ τ τH τR

dimensionless deactivation time [–] reactor path length, free variable, ξ = τ/τR [–] space–time or residence time, free variable [s] characteristic time constant of heat transfer [s] total reaction time [s]

There are many approaches of interest not mentioned here, but an extensive review of reactor network synthesis approaches can be found in Hildebrandt and Biegler (1995). The proposed method, which belongs to the systematic generation methods, can be viewed as a procedure built around the kinetic model. The kinetic model of the phenomena is the crucial basis for this method. The phenomena of importance for the design will include not only the reaction kinetics, but also phenomena like mass and heat transfer, catalyst deactivation, heat generation and consumption, and chemical equilibrium. The structure of the chemical system is defined by the reaction mechanisms taking place, such as parallel and consecutive reactions, but also phase transitions. The generated reactor design will be no better than the predictive power of the kinetic model. Therefore, in order to achieve better reactor designs, better reaction kinetics models are also required, that capture the entire conversion and temperature range. For isothermal conditions, the method is presented in part I of this article (Hillestad, 2004). Here, the method is extended to non-isothermal conditions. A non-isothermal process is a process where the temperature is not necessarily constant along the reactor path and the temperature profile is subject to optimization. There may or may not be heat transfer along the reactor path. This is also called a polytropic process. Since, for most chemical reactors, the temperature is not constant and heat is to be transferred to or from the reaction, a non-isothermal design and operation is the more industrially relevant option. It is possible to let the temperature itself be viewed as a design function. This is, however, a less realistic approach and involves extra constraints to avoid infeasible temperature profiles. For example, a CSTR exhibits constant temperature along the residence time, and unless constraints on the temperature are imposed, the solution will be wrong. Furthermore, in order to achieve the desired temperature profile it is quite likely that the required heat transfer area or coolant temperature will be prohibitive. Therefore, a more realistic

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

approach is by introducing an extra equation representing the steady-state conservation of energy and constructing design functions associated with the heat transfer.

2. Reactor model The exercise of mixing two reaction volumes made in part I of this article can also be made for the non-isothermal case. When two volumes are mixed, the temperature is treated the same way as the concentration vector for the isothermal case, so long as the two volumes have the same specific heat capacity and density.

The temperature equation, constituted by the conservation of energy, also conforms with the dispersion model. The same results derived from the dispersion model for the isothermal case also applies to the non-isothermal case such as the yield and space–time yield. A term representing the heat transfer must be part of the model, and this term is normally proportional to the temperature difference between the reaction and cooling/heating medium: T = T − Tw . This may not always be the case, as when heat radiation is significant. An extra term may be introduced to take radiation into account. Heat transfer through evaporation is included by this formulation, since a chemical component may be modeled in different phases as different species. When a component is transferred from one phase to another, such as A() → A(g), it is essentially equivalent to a chemical reaction with a rate expression and heat of reaction. PFR. For a PFR, the steady-state reactor temperature profile as a function of the space–time, τ ∈ [0, τR ], is: (1)

The ratio A∗ (τ) = (UA)/(UA)R is the relative heat transfer area-coefficient, or in short only area. The parameter U is the overall heat transfer coefficient, being a lumped parameter composed of the heat resistance from an average radial position, through the wall and to the heating/cooling medium. It depends on the film coefficients on both sides of the wall, which again depend on the fluid conditions, in particular the fluid velocity. A∗ (τ) is a relative measure compared to the total heat transfer, (UA)R , which is the sum of all heat transfer areas along the reactor path. A∗ (τ) will always increase (or have zero slope) with the space–time. The initial condition of the differential equation is the feed temperature, TF , or the temperature of the upstream unit. The parameter NSt is the dimensionless Stanton number of heat transfer. The vector h contains the reduced heat of reaction of all reactions defined by the reaction system r: h=

−H cp

(2)

(UA)R q cp

(3)

The Stanton number increases with the reactor size since the total heat transfer area AR is involved. The Stanton number can also be written as the ratio between two time constants: τR NSt = (4) τH The time constant τH is the characteristic time constant of heat transfer: cp τH = (5) (Ua)R a=

2.1. Conservation of energy

dT dA∗ = h r(c, T ) − NSt T dτ dτ

NSt =

1103

AR VR

(6)

The characteristic time constant of heat transfer is independent of the total reactor size. It depends on the specific heat transfer area, a, which is the average heat transfer area over the entire reactor volume. CSTR. For a CSTR, the steady state temperature at variable space–time τ ∈ [0, τR ] is an algebraic equation implicit in temperature: T − TF = τh r(c, T ) − NSt T (A∗ (τ) − A∗0 )

(7)

Here, A∗0 is the relative heat transfer area at the point where the CSTR starts. For a given space–time, τ, the CSTR temperature will be constant over the space–time. But, by letting the space–time be a free variable, the temperature will vary with the space–time according to Eq. (7). The A∗ -function. Considering a single reactor, the heat transfer area will normally increase proportional to the reactor volume, i.e. A∗ = τ/τR . However, when several reactors are encountered along the path, the optimal heat transfer area distribution, A∗ , will likely be a different function. In general, A∗ may be considered to be a smooth function, but that is almost infeasible for practical reasons. A more realistic scenario is to let A∗ be a piece-wise linear function. The slope of each section is determined by the heat transfer area in that section of the reactor path. By definition, A∗ (τ = 0) = 0 and A∗ (τ = τR ) = 1. Fig. 1 illustrates a possible A∗ -function with three different heat transfer sections. For a single CSTR, the A∗ -function will always increase linearly with the variable space–time, because a CSTR will always "see" an average heat transfer area and coefficient. For a single PFR, on the other hand, the A∗ -function may

Fig. 1. A piece-wise linear A∗ -function, where most of the heat transfer area is assigned to the first section, no heat transfer area to the middle section and the rest to the third.

1104

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

be divided into several sections with different slopes or heat transfer area distributions. For the sake of simplicity and as a first approach, the sections of different heat transfer areas are defined to be identical to the path length of a single reactor. That is, each reactor along the path can only have one heat transfer coefficient, UA, increasing proportional to the size of the reactor. This is also the most realistic design target and will be used here in the sequel.

r 

2.2. Component mole balances and conservation of energy

dx dA∗ = R(x) − Ux dτ dτ

(8)

x − xF = τR(x) − Ux(A∗ (τ) − A∗0 )

(9)

U is a diagonal matrix, and in this case it contains only one non-zero element: U = diag(0 · · · 0, NSt ). 2.3. Composite model The algebraic Eq. (9), representing the CSTR model, may be written on differential form. By differentiating the equation with respect to space–time, the CSTR equation can be written:
  
A∗ − A∗0 −1 dx dA∗ R(x) − Ux = I−τ J−U dτ dτ τ (10) As already discussed, a single CSTR will always “see” an average heat transfer coefficient, and therefore the only realistic A∗ is a linear function. The term (A∗ − A∗0 )/τ will therefore always be identical to dA∗ /dτ. Thus, the CSTR equation on differential form may be written slightly different:
 
 dA∗ −1 dA∗ dx = I−τ J−U (11) R(x) − Ux dτ dτ dτ Based on this, two design functions associated with heat transfer are defined. One is the temperature profile of the cooling/heating medium and the other is the relative heat transfer area distribution:

uH =

(12) dA∗

1 τH dξ

dj = 1

(14)

j=1

Defining the state vector x = [c , T ], and a new stoichiometric matrix N augmented with a row h , the augmented component rate vector becomes an extra element, RT = h r [K/s], and is calculated as before, R = Nr. The same formulations as for the isothermal case can now be used. Only an extra term, representing the heat transfer, is added to the last equation. Formulation of the PFR and the CSTR reactor equations, as a function of the space–time τ, are given below:

uT = Tw ,

Then, the term U dA∗ /dτ = EuH , where E = diag(0 · · · 0, 1). The design function uH contains the average heat transfer area and the distribution of area along the reactor path. Since A∗ is a piece-wise linear function, uH is a piece-wise constant function on ξ. Let the piece-wise constant parameterizations of uH be denoted uH,j and the path length of each constant section be dj , then the following must hold:

(13)

r  j=1

dj uH,j =

1 τH

(15)

The number of reactors along the path is r. As already shown in part I of this article, for isothermal operation, the CSTR and the PFR equations may be formulated as one composite equation by introducing the design function, uM , associated with mixing profile. By introducing the dimensionless path length ξ = τ/τR , a model comprising both the CSTR equation (11) and the PFR equation (8) can be constructed: dx = τR [I − τR uM (J − EuH )]−1 [R(x) − ExuH ] dξ

(16)

It is straightforward to verify that when uM is zero, Eq. (16) becomes the PFR equation (8), whereas a value proportional to the path length ξ Eq. (16) becomes the CSTR equation (11). This is illustrated in Fig. 2 by the sawtooth shaped curve. In the same figure, a possible design function uH is also shown. An intermediate slope β of uM , where β ∈ 0, 1, can be interpreted as a CSTR in parallel with a PFR, and both with equal residence times. The volume fractions going through the CSTR and PFR are β and 1 − β, respectively. Distributed feed streams may be added the same way as for the isothermal case, since the temperature of the feed stream is mixed in an similar way as the concentration is mixed. Extra feed streams can only be located where uM = 0, either along a PFR or at the inlet of a CSTR: dx = f = τR [I − τR uM (J − EuH )]−1 [R(x) − ExuH ] dξ + (xF − x)uF

(17)

The composite model (17) contains several design functions to be optimized. The task of putting together a sequence of ideal reactors along the path is represented by the design function uM (ξ), the task of distributing the available feed streams is represented by the design function uF (ξ), and the task of distributing the available heat transfer area among the ideal reactors is represented by the design function uH (ξ). In fact, for each available feed composition xF, , a design function uF, can be assigned, but for the sake of simplicity only one is used here. In addition, a design function uT defining the temperature of the cooling/heating medium is also available, but

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

1105

3.2. Design with catalyst deactivation Most packed bed catalytic reactors suffer from catalyst deactivation. In order to keep the conversion per pass sufficient and to avoid catalyst reloading too often, it is necessary to have some spare catalyst volume with fresh activity. As the catalyst deactivates, the rate of formation profile will be shifted along the reactor path. Let t be the long term time of deactivation, let TL be the catalyst lifetime or time interval between each catalyst reload, and let the dimensionless deactivation time be θ = t/TL . Furthermore, let the original kinetic model containing deactivation kinetics be denoted ˜ x), the extended space–time yield that includes catalyst R(θ, deactivation is:  1 1 ˜ k (θ, x) dθ dξ ϑ= R (20) 0

0

This implies that the average kinetic model over the time horizon between each catalyst loading should be the basis for the design: Fig. 2. Admissible design functions uM and uH as a function of the path length.

 R(x) =

1

˜ x) dθ R(θ,

(21)

0

here it is chosen to be treated as a given profile. The reason is that the cooling/heating temperatures are often determined by the auxiliary steam system and/or available process streams for heat integration.

The rate model to be used for reactor design is the averaged kinetic model over the catalyst lifetime TL . This illustrates how catalyst deactivation can be taken into account in the design by this method.

3. Optimization

3.3. Optimal mixing profile, feed and heat transfer distribution

3.1. Objectives The design criterion applied for the isothermal case may also be applied here. However, for the non-isothermal case an extra degree of freedom is available, cooling or heating. There is a cost associated with heating and cooling and can be compared to the benefit of maximizing the production of the key component. It is easy to modify the criterion to incorporate the cost of heat transfer. However, here the same criterion as for the isothermal case is used. The yield (ψ) and space–time yield (ϑ) of the key component at the end of the reactor path are defined as earlier:  τR ψ = c1,k − cF,k = Rk (x) dτ (18)

With the composite model of Eq. (17) formulated as a system of ordinary differential equations, optimal control theory and the maximum principle may be applied to solve for the design functions uM , uF and uH . Though, other methods may also be applied, such as mixed integer nonlinear programming and advanced dynamic optimization techniques, the maximum principle requires construction of a Hamiltonian, which may be identical to what is chosen for the isothermal case. The space–time yield is chosen to be maximized, subject to having a minimum yield. The space–time yield can be written more conveniently as:  ϑ= 0

0

c1,k − cF,k = ϑ= τR



1 0

Rk (x) dξ

(19)

The space–time yield (ϑ) of a key component is to be maximized because it is equivalent to maximizing the production efficiency. In order to avoid too small conversions, a minimum yield is required. This is explained in more detail in the previous paper.

1

ek f dξ τR

(22)

The vector f = dx/dξ is the composite model defined in Eq. (17), while the vector ek is the kth unit vector. Hence, the Hamiltonian is: H=

ek f + p f τR

(23)

The vector p is the co-state vector or the adjoint vector of the state vector, x. According to the maximum principle, the

1106

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

co-states are given: ∂H dp =− dξ ∂x

the vector uH . Furthermore, the path length of each reactor, denoted dj , are stacked in a vector d: (24)

Since the final states are not included in the objective function, ϑ (Lagrange form), the final value of the co-state vector is zero, p(1) = 0. For the composite model (17), the co-states equation can be derived analytically:   ek dp  = −L p + (25) dξ τR L = [I − τR uM (J − EuH )]−1 
  dJ × τR J − EuH + uM − IuF dξ

(26)

The gradient of H with respect to the design function uM is only slightly different from that of the isothermal case:   dH ek = τR + p [I − τR uM (J − EuH )]−1 (J − EuH ) duM τR (28)

The design function that maximizes the Hamiltonian is the one where dH/duM = 0 at any point along the reactor path. If the gradient along the path is zero, the Hamiltonian will be constant along the path. However, the design function cannot attain any value and the admissible values are the one defined by the CSTR and PFR equations. Therefore, the gradient dH/duM will not be zero and the Hamiltonian will not be constant along the reactor path. The gradient of H with respect to the design function uH can also be derived analytically:   dH ek = −τR + p [I − τR uM (J − EuH )]−1 duH τR × [Ex + uM Ef]

(30)

d = [d1 , . . . , dr ]

(31)

The condition stated in Eq. (15) is on vector notation d uH = 1/τH , and represent necessary conditions for uH to be consistent. The direction in which to go in order optimize the distribution of heat transfer area of reactor j is represented by the average gradient over the reactor: δuH,j =

Eqs. (17) and (25) constitute a two point boundary value problem. The numerical solution of this system can be provided by any two-point boundary value solver, but for the sake of efficiency the fact that the system is block triangular should be exploited. The gradient of H with respect to the design function uF is equivalent to that of the isothermal case:   dH ek = + p (xF − x) (27) duF τR

× [f − (xF − x)uF ]

uH = [uH,1 , . . . , uH,r ]

(29)

Similarly, the design function that maximizes the Hamiltonian is the one where dH/duH = 0 at any point along the reactor path. Also in this case, the design function cannot attain any value. The design function uH is parameterized by a piece-wise constant function, one constant for each reactor along the path. The piece-wise constants, denoted uH,j , are stacked in

 (dH/duH )i i∈Ij

dj

(32)

The vector δuH = [δuH,1 , . . . , δuH,r ] will be used to move uH towards the optimal solution. However, there is no guarantee that this correction leads to a consistent solution defined by condition (15). If the initial value is consistent, the following procedure will guarantee that the updated solution is consistent:  uH := uH + &[δuH − (τH δuH d)uH ]

(33)

The parameter & ≥ 0 is the length of the step, and must be selected such that the heat transfer area in any reactor does not become too small, which may eventually lead to unstable conditions. The optimization algorithm calculates the optimal mixing profile, uM , the optimal feeding profile, uF , and optimal distribution of heat transfer area, uH , at a given total reaction time and overall heat transfer area per volume. Based on the gradients, the Hamiltonian is maximized by using the conjugated gradient method. The constraints on the design function are accounted for by cutting the function so that it becomes an admissible function. The total reaction time (τR ) and the overall heat transfer area (τH ) may be optimized in an outer loop. In the outer loop, the objective may be to maximize the STY, yield or some other economic function. Another issue of importance is to allow the feed temperature to any reactors be optimized. Heat may be exchanged, and feed stream temperatures may be viewed as variables subjected to optimization. This is a matter of available streams for heat integration and a matter of equipment cost. The fourth design function, the temperature of the cooling/heating medium uT , may be optimized in a similar way as the other three design functions. It seems natural to let it be a piece-wise constant function, one temperature for each reactor along the path. However, in the cases below, the function is chosen a given constant and is not optimized. The coolant temperature is often determined by the auxiliary steam system and available streams for heat integration.

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

4. Case studies Two different chemical systems of industrial importance are studied; the methanol synthesis and the steam methane reforming process. In both systems, the bulk phase is gas. Because the number of moles is changing during the reaction and with concentration as state variables, the total pressure will change with the extent of reaction. This is physically incorrect, as the pressure will be equalized due to the momentum balance. Therefore, for gas phase reactions, it is more convenient to use mole fractions instead of concentration as the basis. The concentration vector is the product of the total concentration and the mole fraction vector, c = ctot y. If Rc is the component reaction rate vector, as used previously, with unit [mol m−3 s−1 ], and the new reaction rate vector, R, now is on mole fraction basis with the unit [1/s], and assuming constant total concentration, the conversion is:  (Rc − y Rc,j ) R= (34) ctot The last term is due to the fact that the mole number is changing. With this modification of the reaction rate, c can be substituted with y and all the rest is exactly as described above. The state vector is hence the mole fractions and the temperature, x = [y , T ]. 4.1. Methanol synthesis In the synthesis of methanol from synthesis gas (CO, H2 and CO2 ) on a Cu/ZnO/Al2 O3 catalyst, the three most dominating reactions assumed to take place are: CO + 2H2
CH3 OH

(35)

CO + H2 O
CO2 + H2

(36)

CO2 + 3H2
CH3 OH + H2 O

(37)

The kinetics and thermodynamics of the synthesis is not very favourable in formation of methanol at the required pressure and temperature. Therefore, unconverted reactants are recycled. The effluent stream is cooled and essentially all methanol and water are knocked out in a separator and synthesis gas is recycled. Due to inerts a fraction of the recycle gas will be purged in order to maintain the pressure. A normal thing to do is to let the hot effluent stream be heat integrated with the cold recycle gas before it is fed to the reactor. Here, once through optimization at fixed feed composition is considered. The feed composition, being a mixture of makeup gas and recycled gas, is set to 7 mol% CO, 6 mol% CO2 , 69 mol% H2 , and the rest are inerts. The feed temperature is 210 ◦ C. Thermodynamically it is favourable to operate at low temperature and high pressure, while the kinetics require high temperature. The pressure is set to 82 bar. The kinetic model developed by Graaf, Stamhuis, and Beenackers (1988) is assumed to be true, although more recent research indicate that methanol is formed merely by hy-

1107

drogenation of carbon dioxide. The reader is referred to the article of Graaf et al. (1988) for a description of the kinetics, which will not be repeated here. Here, intra-particle transport resistance is neglected. In addition, catalyst deactivation is not considered. With a deactivation model available including a typical catalyst lifetime, the kinetic model should be averaged as described above. Since it is not available, the applied kinetic model is for a fresh catalyst. Formation of byproducts, such as ethanol and dimethylether, is not part of the kinetic model, but it is known that byproducts formation increase at increased peak temperature. Though, a more accurate study of the reactor design should include the kinetics of byproduct formation and deactivation as well. Case 1. One available methanol synthesis process technology is where a fixed bed catalyst is packed in tubes and where heat is removed through boiling water on the shell side. This is the Lurgi technology described in an article by Supp (1981). The Lurgi design is close to a PFR. This reactor concept has a much enhanced heat removal capacity, because of the large heat transfer area. Based on information of a typical temperature profile, the characteristic time constant for heat removal is approximately τH = 0.7 s and the total reaction time is set to τR = 14 s. By keeping τR and τH constant at these values and the coolant water temperature constant at 250 ◦ C, the optimal fluid mixing profile is found to be a single PFR. It is also found that it is not beneficial to distribute parts of the feed along the PFR. This confirms that the Lurgi design is close to being optimal with respect to fluid mixing, given the large heat transfer area. Fig. 3 shows the concentration and temperature profile at these conditions. Case 2. Heat transfer area is expensive, so it is of interest to study the consequences of reducing the total heat transfer area. What will be the optimal mixing profile when the heat transfer area is reduced by a factor of, say, 3? The characteristic time constant for heat transfer is then tripled, τH = 2.1 s. The total reaction time is kept at 14 s and the cooling water temperature is 250 ◦ C as in Case 1. In this case, the heat transfer area is equally distributed over the reactor path. The peak temperature will naturally increase at reduced heat removal capacity. Increased temperature is thermodynamically unfavourable, as the synthesis reactions are exothermic. The optimal mixing configuration turns out to be a PFR–CSTR– PFR sequence as shown in Fig. 4. With a CSTR in the middle, the peak temperature is reduced. Fig. 5 shows the concentration and temperature profile with this mixing configuration. Still the peak temperature is almost 10 ◦ C higher than in Case 1. The methanol yield is 6.51 mol%, whereas 6.71 mol% in Case 1. Also here, it is found that it is not beneficial to distribute any feed along the reactor path. Case 3. This case is with reduced heat transfer area as in Case 2, but where the heat transfer area distribution is optimized. The total reaction time is still 14 s and the cooling water temperature is 250 ◦ C as in Cases 1 and 2. The optimal mixing configuration turns out to be a CSTR–PFR sequence as shown in Fig. 6. The residence times of the CSTR and the PFR are 4.06 and 9.94 s, respectively. Fig. 7

1108

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

Fig. 3. Case 1: the methanol synthesis solution at high heat removal capacity.

shows the concentration and temperature profile along the path with the CSTR–PFR sequence, and Fig. 8 shows the heat transfer distribution of the two reactors. The optimized heat transfer area density is less in the CSTR, and the corre-

sponding characteristic time constants of heat transfer are 5.58 and 1.67 s for the CSTR and the PFR, respectively. The same again, it is not beneficial to bypass any feed to the PFR.

Fig. 4. Case 2: the optimal mixing profile at reduced heat removal capacity.

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

1109

Fig. 5. Case 2: the optimal solution at reduced heat removal capacity.

The peak temperature is approximately 7 ◦ C higher, while the methanol yield is slightly higher than that of Case 1. Provided that the increased peak temperature does not lead to too much byproducts formed, this may be a very good design

alternative to that in Case 1. The CSTR may be realized by a fluidized bed reactor and the PFR by a packed bed tube reactor. Assuming the benefits of reduced total heat transfer area is not eaten up by the cost of having two reactors instead

Fig. 6. Case 3: the optimal mixing profile at reduced heat removal capacity and optimal heat area distribution.

1110

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

Fig. 7. Case 3: the methanol synthesis solution at reduced heat removal capacity and optimal heat area distribution.

of one, this alternative may also be more cost effective. It is anticipated the benefits are more pronounced for large-scale plants. The practical issues of introducing a fluidized bed reactor in the synthesis loop is not considered here. With this structure, it is possible to reduce the peak temperature by cooling the intermediate stream or by using a lower coolant temperature in the first reactor. Depending on the available pressure levels of the auxiliary steam system or available streams for heat integration, different coolant temperatures of the two reactors may be possible. 4.2. Steam methane reforming An important process for production of synthesis gas or hydrogen is the primary steam reformer. Typical data of

the process and a proposed kinetic model are described by Froment and Bischoff (1990). Methane and steam are fed at 520 ◦ C and 29 bar. Due to recycling, the feed composition is 21.23 mol% CH4 , 71.30 mol% H2 O, 0.21 mol% CO, 2.59 mol% H2 , 1.19 mol% CO2 , and 3.48 mol% N2 . The three dominating reactions assumed to take place are: CH4 + H2 O
CO + 3H2

(38)

CO + H2 O
CO2 + H2

(39)

CH4 + 2H2 O
CO2 + 4H2

(40)

The first and the third reactions are strongly endothermic so the conversion of methane to hydrogen requires much heat. The kinetic and thermodynamic model can be found in Xu

Fig. 8. Case 3: the optimal heat transfer area distribution at reduced heat removal capacity.

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

1111

Fig. 9. Case 4: the concentration and temperature profile of a PFR steam reforming tube.

and Froment (1989), and is not repeated here. An approximate catalyst pellet effectiveness factor of 0.03 is used for all three reactions. The key product component is hydrogen. Case 4. A conventional steam reformer described in Froment and Bischoff (1990) have many tubes in parallel filled with a Ni-type catalyst. They are mounted in a furnace in order to transfer the required heat. The characteristic time constant for heat transfer is set to τH = 1.5 s, which is quite low. The reaction time is set to τR = 4 s. The furnace temperature is set to 825 ◦ C and assumed constant along the path. The optimal fluid mixing configuration of the stream reforming process is a PFR. Distributed feed is of no benefit, which is more or less evident in this case. A conventional steam reformer tube is close to a PFR, which is in accordance with the result of calculations made here. Fig. 9 shows the concentration and temperature profile along a PFR steam reformer tube. By doubling or by halving the heat transfer area, the optimal mixing profile is still a PFR.

5. Discussion and conclusions A new method for generation of conceptual reactor designs is proposed. It can be considered as an alternative to the existing methods such as the attainable region method and superstructure optimization. The method is based on a new composite model formulation comprising both the CSTR and PFR models. A design function uM defining the sequence and

size of ideal reactors and another defining the distributed feed along a path, uF , are targeted. In this paper, non-isothermal conditions are treated and the model formulation is augmented with the temperature equation. This gives an extra design function, uH , the distribution of heat transfer. The design functions are solved by using optimal control theory and the maximum principle. The composite model formulation and how it is parameterized by the design functions are new concepts. It constitutes a concise and realistic way of formulating the reactor design problem. Though, the main focus here has been on presenting the method, it facilitates the focus on the phenomena taking place. The design engineer may concentrate on developing a model of the rates at which the phenomena proceed, and specifying available feed compositions and temperatures. The method, implemented as suitable computer software, will automatically generate the design functions. These design functions are translated to a reactor flowsheet. Two industrially important processes are studied; the methanol synthesis and the steam methane reforming process. It is demonstrated that the heat transfer area can be reduced substantially in a methanol synthesis reactor without reducing the methanol yield, by introducing a CSTR–PFR sequence and optimal distribution of heat transfer area. It is also found that for the steam methane reforming system, the optimal mixing profile is a PFR. In addition to the items mentioned in part I of this article, further development of the method should include the pres-

1112

M. Hillestad / Computers and Chemical Engineering 29 (2005) 1101–1112

sure as a design function. A natural choice is to let the pressure be a piece-wise constant function, one pressure level for each reactor. The fact that there may be a pressure drop across a reactor can be taken into account, but may also be neglected as a first approach. The pressure as a design function is particularly interesting for systems where phase transitions occur but also for pure gas systems. The proposed method may then be applied to the design of reactive distillation systems. Furthermore, dilution of the catalyst concentration may also be an appropriate design function in some systems as described by Hwang and Smith (2004). For other chemical systems it might be beneficial to use different catalysts along the path, as described by Wheeler, Jhalani, Klein, Tummala, and Schmidt (2004). The selection of catalysts may then be represented by one or several design functions. Still more, membrane reactors, where the mass transfer flux is proportional to the concentration difference, is readily adapted to this method in a similar way as the heat transfer. The design function may be the distribution of membrane area. References Balakrishna, S., & Biegler, L. T. (1992). Constructive targeting approaches for the synthesis of chemical reactor networks. Industrial and Engineering Chemistry Research, 31, 300–312.

Feinberg, M. (2002). Toward a theory of process synthesis. Industrial and Engineering Chemistry Research, 41, 3751–3761. Froment, G. F., & Bischoff, K. B. (1990). Chemical Reactor Analysis and Design. New York: Wiley. Graaf, G. H., Stamhuis, E. J., & Beenackers, A. A. (1988). Kinetics of low-pressure methanol synthesis. Chemical Engineering Science, 43, 3185–3195. Hildebrandt, D., & Biegler, L. T. (1995). Synthesis of reactor networks. AIChE Symposium Series, 91, 52. Hillestad, M. (2004). A systematic generation of reactor designs-I. Isothermal conditions. Computers & Chemical Engineering, 28, 2717– 2726. Hopley, F., Glasser, D., & Hildebrandt, D. (1996). Optimal reactor structures for exothermic reversible reactions with complex kinetics. Chemical Engineering Science, 51, 2399–2407. Hwang, S., & Smith, R. (2004). Heterogeneous catalytic reactor design with optimum temperature profile. I. Application of catalyst dilution and side-stream distribution. Chemical Engineering Science, 59, 4229– 4243. Jackson, R. (1968). Optimization of chemical reactors with respect to flow configuration. Journal of Optimization Theory and Applications, 2, 240– 259. Supp, E. (1981). Improved methanol process. Hydrocabon Processing, 3, 71. Wheeler, J., Jhalani, A., Klein, E. J., Tummala, S., & Schmidt, L. D. (2004). The water-gas-shift reaction at short contact times. Journal of Catalysis, 223, 191–199. Xu, J., & Froment, G. F. (1989). Methane steam reforming, methanation and water-gas shift. I. Intrinsic kinetics. AIChE Journal, 35, 88– 96.