distillation processes using detailed dynamic models

Computers them. Engng

Pergamon

0098.1354(95)00170-0

Vol. 19, Suppl., pp. S167-S174.1995 Copyright 8 1995 Elsevier Science Ltd Printed in Great Britain. All rights reserved 0098-1354/95 $9.50 + 0.00

SYNTHESIS OF BATCH REACTION/DISTILLATION PROCESSES USING DETAILED DYNAMIC MODELS

MS. CHARALAMBIDES, N. SHAH AND C.C. PANTBLIDBS~ Centre for Process Systems Engineering, Imperial College, London SW7 2BY, UK.

ABSTRACT This paper considers the synthesis of processes comprising batch reaction and distillation tasks. Complex process features, such as material recycles and shared intermediates are taken into account. The problem is formulated as a multistage optimal control problem. A control vector parametrisation approach is utilised to reduce this infinite dimensional dynamic problem to a finite dimensional nonlinear programming problem. An example involving an exothermic reaction step followed by a batch distillation is presented to illustrate the proposed approach.

KEYWORDS

Batch processes; process synthesis; batch process optimisation; optimal control.

INTRODUCTION Batch processes implement recipes involving sets of physical, chemical and biological transformations (“processing tasks”) that transform the raw materials to the final products, usually via intermediate materials. In recent years, a substantial body of work has been reported in the literature on the problems of scheduling the operation of existing batch plants, and designing new ones (see, e.g., the reviews by Reklaitis 1989, 1992). In this work, it is usually assumed that the process recipe is fixed and known a prioti. This information includes both the structure of the processing network and the precise operating policies for the individual processing tasks, as well as the properties of all materials involved in the process. Batch process synthesis is concerned with the design of the process itself using appropriate models for the individual processing tasks. In the most general case, one seeks to establish both the set of processing tasks that will transform a given set of raw materials into a given set of final products, and the operating policy for each task. This paper is concerned with a more restricted form of the problem which assumes that the structure of the processing network is known and seeks to determine appropriate operating policies for the various tasks in it. This is by no means a trivial problem: most batch unit operations have one or more control variables (e.g. feed addition rate, reflux ratio) that can be manipulated, and the operating policy involves not only the appropriate batchsize and processing time for each task, but also the time variation of the available controls. The interactions between process synthesis and plant design for a fixed task network were considered by Iribarren (1985) who employed simplified algebraic unit operation models to determine optimal operating conditions at the preliminary design stage. A decomposition approach to this problem was proposed by Barrera and Evans (1989), and this was extended to more general processes by Allgor et al. (1994). The interactions between unit operation and plant scheduling were considered by Salomone and Iribarren (1992) and by Tricoire (1992) who demonstrated the advantages of a simultaneous approach to the two problems over decoupled and iterative techniques. Again, simple algebraic models relating the unit operation performance to t

Author to whom correspondence should be addressed. S167

S168

i European Symposium on Computer Aided Process Engineering-5

key decision variables (e.g. batchsizes) were employed. A major deficiency of these approaches is that either they rely on tamer simple models for describing what are essentially complex time-dependent unit operations, or they utilise a decomposition approach to the problem, where detailed models may be used at some level of the hierarchy but the relationship between optimisation objectives at different levels and the overall process objective is not always clear, This paper advocates a simultaneous approach to the problem of determining optimal operating strategies for entire batch processes, and illustrates its application to a reaction/separation system. In contrast to previous work, we employ detailed dynamic models for the individual unit operations, thus ensuring the feasibility and optimality of the solutions obtained within the modelling error. The basis for the work reported here is the batch process synthesis formulation presented by Charalambides et al. (1993). This is reviewed in the next section. The subsequent section &fines the process synthesis problem being addressed here. This is followed by a detailed example involving batch reaction and distillation operations with multiple recycles of material. The paper concludes with some remarks concerning the applicability of the work presented.

CHARACTBRISATION OF BATCH PROCESS RECIPES We assume that the structure of the process being designed is given in the form of a State-Task Network (SIN) (Kondili et. al, 1993). An example of such a STN is shown in Fig. 1 which represents a process involving a reaction and a distillation operation. The latter produces four different products, two of which are recycled. In the standard STN notation, circles represent states, each corresponding to a different type of material while rectangles represent &z&s,each transforming a set of input states to a set of output states.

-I Fig. 1. State-task network representation of process recipe Each state s in the SIN is character&d by a set of intensive properties zs. Some of these may be given (e.g. those corresponding to the raw material states), but in general these are variables to be determined by the process synthesis algorithm. Typically, commercial and operability considerations may dictate that these lie within given bounds: Z~IZ15Z~

vs

(1)

The operation of most types of batch task is described by mixed systems of differential algebraic equations (DABS) of the form: fi(xi(r)* ki(r), vi(r), a,(r). vi) = 0

Vi. tG[Ov Sl

(2)

where Xi(t) are the differential variables, yi(t) the algebraic variables, a,(t) possible control manipulations and Vi time invariant parameters (e.g. equipment sixes) respectively. Both the time-invariant parameters and controls are likely to be bounded from above and below. Additionally, safety, operability and other considerations may result in path constraints that must be enforced throughout the operation of the task, and point constraints that must hold at a finite set of time points. In many real applications, the duration +i of each task i is determined implicitly through end-point constraints. Each task i consumes a set of states Si and produces another set of states gi. The initial condition of the task is determined by the intensive properties of its various input states and the precise amount of each state used by the task These amounts are also variables to be determined by the optimisation. Similarly, the amount and properties of each state produced by a task will usually depend on the final condition of the task.

European Symposium on Computer Aided Process Engineering-S

S169

PRELIMINARY PROCESS SYNTHESIS ALGORITHh4 The synthesis problem considered here involves the design of a new plant to implement a process of a given

STN structure. We seek to establish the characteristics of the recipe, the size of the equipment items utilised, and the operating policy of each task to meet the required product specifications. It is assumed that unlimited intermediate storage is available for each state. Furthemore, all batches of the same task take place in a single dedicated equipment item and are produced in exactly the same way over the time horizon of interest. A complete mathematical statement of the problem is given by Charalambides et al. (1993). The optimisation variables are the intensive properties of each state, the processing time, control variables and time-invariant parameters for each task, and the amounts of each input state used by each task. Because of the time variation of the control variables, in general this is an infinite dimensional dynamic optimisation problem. A control vector parametrisation algorithm (Vassiliadis era!., 1994a) is used to reduce it to finite dimensional form which can be solved using standard nonlinear programming (NLP) algorithms. Objective function and constraint gradients are calculated by integrating the system sensitivity equations, and inequality path constraints are handled using the hybrid technique proposed by Vassiliadis ef al. (1994b).

EXAMPLE The process under consideration is described by the STN of Fig. 1 and is to be implemented in a batch plant of the form shown in Fig. 2. Three tasks are utilised to transform pure feedstocks, Feed4 and FeedB to the final product, Pr. Three additional materials are produced. One of them, Res, is too heavy to be of any practical use within the process, but the other two can be recycled. Obviously, when the plant is being started up, the recycles are not yet avaiIable. In fact, even when they do become available, a number of batches may have to be processed before the plant reaches a periodic steady state. However, since we are interested in optimising the long term operation of the plant rather than its start up phase, we consider this periodic state directly by assuming that the material recycled into the reaction and the mixing/heating steps is of exactly the same quality as that produced by the distillation.

-r

--__-_______I

Fig. 2. Plant flowsheet Process Description

The main reaction A + B + C and the side reaction B + C + D take place in the liquid phase in the reactor. The side reaction is favoured at higher temperatures and consumes the desired product C. The reactor contents must be kept below a maximum temperature of 350 K in order to avoid vapourisation of the reacting mixture; temperature may be controlled by manipulating the flowrate of the cooling water in a cooling jacket.

s170

European Symposium on Computer Aided Process Engineering-5

The reactor operation is modelled in terms of a set of standard material and energy balances. The cooling driving

force is taken as the difference between the outlet cooling water temperature and the temperature of the reactor contents which are assumed to be well-stirred. The maximum cooling water flowrate is 3.6 tonslhr, but the actual flowrate can be varied in a piecewise constant fashion over the duration of the reaction. For the purposes of the optimisation, we assume that this is done over 5 time intervals of variable length. The mixing-heating task pre-heats the reactor product Intl to 360 K, a temperature close to the bubble point of the material fed to the distillation column. At the same time, it mixes Intl with any amount recycled from state Red. The heating medium used is saturated steam at a pressure of 4 bar. The steam flowrate is assumed to be constant throughout the heating operation, up to a maximum of 1 tonlhr. The distillation column comprises 12 trays and is used to produce three distillate cuts. Saturated steam at 4 bar is used in the reboiler. The volatility of the four components in the process decreases in the order A > B > C > D. Consequently, the first cut Reel recovers the unreacted reactants A and B; clearly these can be recycled back to the reaction. The second cut Rec2 comprises a mixture of the reactants with the desired product C; this may be mixed with fresh reactor product to undergo distillation. The desired product C is recovered in the third cut with a specified minimum purity of 95 %. All three distillate cuts are cooled down to 300 K, using cooling water at 288 K. The residue Res remaining in the column (including reboiler, trays and condenser) at the end of the distillation is primarily component D and is discarded. The operation of the distillation column is modelled assuming constant liquid hold-ups on the trays and the condenser, negligible vapour hold-up, and fast energy dynamics. The liquid holdup in the trays and condenser is 0.05 kmol, and the initial charge to the column fills the reboiler, trays and total condenser with liquid at the charge composition. The separation takes place at atmospheric pressure assuming an ideal mixture. The vapour pressure of each component is calculated through an Antoine equation. Both reflux ratio and reboiler load are manipulated in a piecewise constant fashion. The reflux ratio is allowed to vary between 0.5 to 0.99 and the reboiler load between 550.0 and 1250 MJlhr. For the purposes of the optimisation, the operating time of distillation is divided into 6 intervals of variable duration, with the first and second cuts being produced at the ends of intervals 2 and 4 respectively. The physical data utilised in the reaction, mixing and distillation tasks are provided in the Appendix. Available Equipment It is assumed that the plant will comprise one equipment item dedicated to each task. Equipment capacity is constrained to lie within the limits shown in the first two columns of Table 1. Table 1. Plant equipment characteristics.

Equipment Item Reactor Mixer-Heater Distillation Col.

Minimum Capacity (kmol) 15.0 15.0 15.0

Maximum Capacity (kmol) 50.0 50.0 40.0

Cost Coefficient (rculkmol”.8. hr) 0.086 0.033 0.185

The unit values of the various materials in the process are shown in Table 2l. Note that the values of Int 1, fnt2, Ret 1, Rec2 and Res are negative, reflecting tbe undesirability of building up stocks of these materials. The unit costs for the cooling water and steam utilities are 0.2 and 10.0 rculfon respectively. For simplicity, we assume that the annualised capital cost of all equipment items is given by ai Butchsiz.e”.8. The values of the coefficients ai are given in the last column of Table 1. The objective function to be maximised is the nett revenue per unit time of operation. It comprises the value of the products, the cost of the raw materials and utilities, and the (annualised) capital cost of the equipment.

’ 1124 = Relative CostUnit

European Symposium

on Computer Aided Process Engineering-S

s171

Table 2. Unit values for process materials

Material

Unit value (rculkmol)

FeedA FeedB Intl Int2 Ret 1 Rec2 Res Pr

3.0 3.0 -0.1

-0.1 -0.1 -0.1 -0.1 30.0

Results and Interpretation

The overall model of the process involves a total of 420 differential and algebraic equations describing the task operation. The number of decision variables manipulated by the optimisation is 76. An economic breakdown of the solution obtained is provided in Table 3. It can be seen that the intermediates Intl, Int2, Reel and Rec2 are not allowed to accumulate, but some accumulation of the waste Res is unavoidable. The characteristics of the recipe are summarked in the annotated ST’N of Fig. 3. The quantity (in kmol) shown on each arrow entering a task rectangle denotes the amount of the corresponding state used to make up a batch of this task. The quantity on each arrow leaving a task rectangle denotes the amount of the corresponding state produced in every batch. The actual processing time for each task is shown as r within the task rectangle. For the reaction and mixingiheating tasks, some idle time is introduced between successive batches in order to ensure that the processing rates of all tasks balance correctly. Table 3. Economic breakdown of optimal solution Equipment

Capacity

cost

Item

(km00

(rculhr)

Reactor 39.08 Mixer-Heater 35.71 Distillation Column 32.29 Total equipment cost Material Nett productionrate FeedA FeedB In11 Int2

(kmollhr) -5.313 -7.285

1.61 0.58 2.98 5.17

Revenue (rculhr) -15.94 -21.86

0.000

0.00

0.000

0.00

Reel

0.000

0.00

Rec2 RtT Pr

0.000

0.00

2.284 5.151

-0.23 154.53

Tot01production revenue Utility Consumptionrate (tonlhr)

34.72 Cooling water Steam 0.609 Total utility cost TOTAL NET REVENUE

116.50 cost (rculhr)

6.94 6.09 13.03

98.30

The profile of the cooling water flowrate in the reactor and the resulting temperature profile are shown in Figs. 4 and 5 respectively. Also Fig. 6 shows the optimal variation of the normalised reflux ratio’ in the distillation column. At the optimal solution, the reboiler heat load is at its maximum value of 1250 MJlhr throughout the distillation. The resulting composition trajectories of the various distillate cuts collected in the accumulator are shown on Fig. 73. It is worth noting that the amount of off-cut produced and recycled to the mixer/heater is held to a very low value: it is more profitable for more material to be included in the first main cut which is recycled 2 Defined as the ratio of the flowrate of the reflux stream to that of the vapour stream entering the condenser.

S172

European Symposium

on Computer Aided Process Engineering-5

________. : xA.l.om: :

i

:x

_______ _.

: x* .0.40: j xs .o.m:

x,.o.od

xc.o.om~

0.0!XQ: I T.m.Olt: ‘-‘I----’

’

:xc

-0.255:

:xD

-o.lm:

+OO.OK: -A----*

:%I rO.lfl :T.ZCO.?K: ._______-

________.

i. 1 ._______s

Fig. 3. Optimal process recipe to the reactor.

Mm

054

I.,“,

nnu OR)

IS0

2m

0.m

050

Time

Fig. 4. Optimal cooling water flowrate variation (todhr)

Irn

1.54

21”)

ou,

Fig. 5. Reactor temperature (K)

“3 0.8

mmtMaIn cti

0.7 0.6 0.5 0.4 0.3 02 0.1

: “....

0.0

Fig. 6. Optimal reflux ratio variation

,...’ ,..” _.,.’

xc

,, ..(__..,....

./

*.

Fig. 7. Accumulator composition

trajectories

’ The non-monotonic bchaviour of these variables near the start of each cut is not significant. It is simply due to the arbitaryvalues used to initialise the composition of the empty accumulator.

European Symposium on Computer Aided Process Engineering-S

s173

CONCLUDING REMARKS This paper has demonstrated the capabilities of the batch process synthesis approach originally proposed by Charalambides et al. (1993). Specific emphasis was focussed on processes involving reaction and distillation steps with multiple recycles of material. However, the synthesis formulation is entirely general and can be applied to any process for which dynamic task models are available. The work presented here can be viewed as a generalisation and extension of earlier work on the optimisation of individual batch unit operations (e.g. batch reaction; batch distillation) that relies on the use of dynamic optimisation techniques. Of particularly note is the general and flexible manner in which recycles of material are handled. Batch process synthesis algorithms based on the use of algebraic models, such as that by Salomone and Iribarren (1992) and Tricoire (1992), are a special case of the approach presented here, with equations (2) being replaced by the simpler form: fi(yi,Vi) = 0

Vi

(23

where the processing time ri is a given function of the time invariant parameters Vi. In fact, our implementation allows the mixing of differential/algebraic models for some tasks with purely algebraic models for others. Finally, it is worth pointing out that, because of the use of local optimisation algorithms for the solution of the NLP, the global optimal@ of any solution obtained with our approach cannot normally be guaranteed. This is a common deficiency of optimisation-based design methods which can only be overcome by the adoption of global optimisation techniques.

REFERENCES Allgor, R.J, M.D. Barrera, PI. Barton and L.B. Evans (1994). Optimal Batch Process Development. PIVC. Fifrh Int. Symp. on Process Systems Eng., 153-157, Kyongju, Korea. Barrera, M.D. and L.B. Evans (1989). Optimal Design and Operation of Batch Processes. Chem. Engng. Commun. 82,45-f%. Charalambides, MS., N. Shah and C.C. Pantelides (1993). Optimal Batch Process Synthesis. Paper presented at AIChE Annual Meeting, St. Louis, Missouri. Iribarren, 0.A (1985). Butch Chemical Processes Design. PhD Thesis, University of Massachussetts. Kondili, E., C.C. Pantelides and R.W.H. Sargent (1993). A General Algorithm for Short-term Scheduling of Batch Operations. Part I - MILP Formulation, Comp. them. Engng., 17,21X-217. Reklaitis, G.V. (1989). Progress and Issues in Computer Aided Batch Process Design. Proc. Third Intl. Conf on Founuiztions of Computer-Aided Process Design, Snowmass, Colorado, 214-276. Reklaitis, G.V. (1992). Overview on Scheduling and Planning of Process Operations. Prvc. NATO ASI on Butch Proc. Syst. Engng., Antalya, Turkey. Salomone, H.E. and 0-A. Iribarren (1992). Posynomial Modeling of Batch Plants: A Procedure to Include Process Decision Variables. Comp. them. Engng. 16,173-184. Tricoire, B., (1992). Design and Scheduling of Multiproduct Batch Plants with Application to Polymer Production. PhD Thesis, University of Massachussetts. Vassiliadis, VS., R.W. Sargent and C.C. Pantelides (1994a). Solution of a Class of Multistage Dynamic Optimisation Problems. 1. Problems without Path Constraints. Ind. Eng. Chem. Res. 33, 2111-2122. Vassiliadis, VS., R.W. Sargent and C.C. Pantelides (1994b). Solution of a Class of Multistage Dynamic Optimisation Problems. 2. Problems with Path Constraints. Ind. Eng. Chem. Res. 33,2123-2133.

APPENDM Reaction Kinetics

muin reaction

rate = kl CA Ce

(kmollhrlm”)

side reaction

rate = kz Ce Cc

(kmollhrlm3)

European Symposium

s174

on Computer Aided Process Engineering-5

where CA, CB, Cc and CD are the concentrations of components A, B, C and D respectively. The kinetic rate constants (m”lhr/kmol) are of the Arrhenius form: -17cOo k,

=lOO.OeRIT

(A.l)

-2oGuo

k2 = 120.0 e ‘GT

(A.2)

where T is the reaction temperature (K) and R, the universal gas constant (8.314 kJ/kmoflK). reaction for the main and side reactions are -60000.0 and -50000.0 kJlkmo1 respectively.

The enthalpy of

Cooling Jacket The cooling jacket in the reactor vessel has a surface area of 30 m2 and an overall heat transfer coefficient of 5088.16 kJlm21Klhr is assumed for the heat transfer between the reaction mixture and the cooling water. The inlet temperature of the cooling water is 288 K with a heat capacity of 4200.0 W/ton/K. Enthalpy and Density

Ideal mixture behaviour is assumed throughout. The liquid and vapour specific enthalpies (Htiq and H,, kJ/kmol) of pure components are given in the form: Hli, = a,(T H w,, = a,(T

where AH,

in

- 296)

- 296) + AH,

(A.4)

is the enthalpy of vapourisation.

Table A. 1. Enthalpy Data

Component a, (Wlkmol K) AH,,(kJlkmot)

A

B

C

D

172.3

200.0

160.0

155.0

31000.0

26000.0

28000.0

34000.0

The density of components A, B, C and D in the liquid phase are 11.0, 16.0, 10.4 and 10.0 kmollm” respectively. Vapour Pressure

The Antoine equation for ideal gases is utilised to calculated the vapour pressure (Pw in bar) at a given temperature (T): InPw=-t+b

64.5)

Table A.2. Vapour Pressure Data

Component a (K) b

A

B

C

D

4142.7500 11.7158

3474.5600 9.9404

3741.8800 9.9915

4543.7 100 11.2599