- Email: [email protected]
Dynamic simulation of multicomponent batch rectification with chemical reactions
Computers & Chemical Engineering, Vol. IO, No. 4, pp. 389-398, 1986 Printed in Great Britain. All nghts reserved
Copynght 0
0098-l 354/86 $3.00 + 0.00 1986 Pergamon Journals Ltd
DYNAMIC SIMULATION OF MULTICOMPONENT BATCH RECTIFICATION WITH CHEMICAL REACTIONS P. E. CUILLE Institut du Genie
Chimique,
Institut
National
Polytechnique
de Toulouse,
Toulouse,
France
G. V. REKLAITIS School
of Chemical
Engineering,
Purdue
University.
West Lafayette,
IN 47907,
U.S.A.
(Received 6 February 1985; revision received 29 October 1985; received,for publication 20 December 1985) Abstract-A model is formulated and solution strategies are investigated for the simulation of a staged batch distillation unit with chemical reactions in the liquid phase. The implementation of the model involves three computation phases: calculation of column profiles at total reflux, calculation of initial derivatives of the algebraic variables and integration of the dynamic model itself. The latter calculation is performed by solving the system of algebraic and differential equations simultaneously using an implicit integration method. The possibilities of vectorizing the model are considered and results are presented on the use of parallel computation. Examples with and without reactions on the plates are reported. Scope-For centuries, batch processes were used in a very wide field of application, from food processing to coal and petroleum engineering. Beginning at about World War II. as it became necessary to cheaply and quickly produce large quantities of products, continuous processes became the standard. More recently, because of the development of new processes for the production of small amounts of products with high added value, batch processes returned to fashion. Publications about the mathematical modeling of batch distillation operations without chemical reactions are found in the literature, as early as the sixties. One of the first rigorous mathematical models was proposed by Meadow [l]. The solution of such a model, using various numerical methods, has been presented by Distefano [2]. During the same period Gear [3] published his first work on a numerical method for the integration of ordinary differential equations. More recently, because of the increasing interest in batch processes, some new works were published: Boston et al. [4] proposed the program “BATCHFRAC” as a general simulation package for batch rectifications. This program uses the model proposed by Distefano, with some numerical approximations, and employs the integration routine “EPISODE” [5], based on the numerical method proposed by Gear. The key approximation employed in Boston’s work consists of the calculation of the time derivative of the enthalpies via a low-order difference approximation. Gallun and Holland [6] proposed a more rigorous approach, using contemporary improvements in the numerical solution of sets of coupled algebraic and differential equations, resulting from Gear’s works. Gallun’s approach is thus based on the direct solution of the balance equations (differential) and of the constraint equations (algebraic) describing the system, He also uses a hydrodynamic model for the calculation of the pressure drop between each plate and of the volumetric holdup on each plate. The plates are assumed to be ideal stages for the vapor-liquid equilibrium calculation. Reactions on the plates are not included in this model. In the present case we do consider reactions and nonideal stages, and also examine the possibilities for vectorization. However, neither this work nor the previously cited references consider the dynamics of the column heat-up or start-up phases which result in the attainment of the total reflux condition. Conclusions and Significance-The simultaneous solution of the model algebraic and differential equations by means of the Implicit integrator LSODI proves to be a highly effective and reliable approach for batch rectification with reaction. Calculation of the initial, total reflux profile can be time consuming if the initial reboiler charge is not large compared to the holdup on the plates. However, a simple physically meaningful relaxation approach provides satisfactory convergence in all test cases examined in this work. The problem of estimating the initial derivatives of the algebraic variables could be very reliably resolved by introducing a staightforward polynomial transition function. In general, the computation time required to simulate a column of a given number of stages over a given simulation time interval increases significantly but manageably with number of components and number of reactions. The parallel computation trials showed that the principal benefits of this mode of computation are achieved if the linear equation solver used in the corrector step is recoded in vectorized form. With such modifications, some 40% can be saved in execution time. Use of automatic vectorization yielded only minor execution time savings (3%). Further significant gains due to vectorization can only be expected if properties evaluations can be carried out in vector-compatible form.
THEORETICAL
MODEL
In this section we consider the mass- and energybalance equations for the batch rectification column
consists of the total mass balance and mass balance on each component. The energy-balance differential equation merely consists of the applica-
shown
tion
in Fig.
1. The
set of mass-balance
equations
differential 389
of
the
first
law
of
thermodynamics.
The
390
P. E. CUILLE and G.
V.
REKLAITIS
and dA,.H,, L dt
f-4
(3)
Plate j, 1=V,+,+L,-,-V,-L,+An;A,, dt
A0
*
I
dJ&= dt
w
0
‘0
l_
= V,.H,,,-(L,+D).H,,~+Q,;
(4)
V,+l'Y,+l+Lj-,'X,~,-V,Y, -
L;X,+R;A,
(5)
and w=
V/+I
dt
i
0,
. ffv., + I +Ljm,‘ffl,,m,
- v,.Hv,,--,.4.,$-Q,;
(6)
Reboiler dA, = L,-dt dA,.X, ~ dt
T/,+An,.A,,
=L;X,-
(7)
V,,.Y,+R,.A,
(8)
and dA,. H1.b ~ = L;H,,” dt Fig. 1. Column schematic.
equations following
presented hereafter assumptions:
are valid
under
the
1. The vapor-phase holdup is assumed to be negligible compared to the liquid-phase holdup on each plate. 2. Chemical reactions in the vapor phase will be neglected. 3. The initial state of the column is the steadystate total reflux condition with no reactions. 4. The liquid volumetric holdups on the plates will be assumed to be constant. Thus, our model is directed at simulating the dynamics of the main production period during which the hydrodynamic conditions are not widely varying. 5. The pressure drops and the plate efficiencies are constant during the operation. The following
differential
balance
equations
result.
Condenser dA, -=V,-L,--D+An,.A, dt = V,An,.A,--D(R dA,.X, -= dt
V,.Y,-L,.X,-D.X,+R,.A,
(1) (2)
(9)
+ Q,,.
An analysis of the degree of freedom of this system readily shows that with these differential equations alone the problem is underspecified. In order to obtain a fully determined system, the variables appearing in these balance equations must also satisfy the following equations. Constraint
on the volume
The volume of the liquid phase in the condenser and on each plate is assumed to be constant (we do not use a hydrodynamic model). Consequently the corresponding molar holdups are functions only of the temperature, the pressure and the composition: A, = r, P,IM,
and
P/ = P CX,, T,, f’,).
(10)
This holdup must be equal to that calculated by the integration of the total mass-balance equations. It should be noted in passing that the inclusion of tray hydrodynamic models such as the Francis-Weir equation merely results in additional equations which can be appended at this junction [e.g. 61. The numerical solution approach discussed in the next section will not be materially altered. Constraint
+ l),
- I”,.&,
due to vapor-liquid
equilibria
For given values of the temperature, and the composition: for each plate j, Y,=Y,+,.(l
-e,)+e;K;X,
the pressure
(11)
391
Simulation of a staged batch distillation unit system of linked differential and algebraic has to be solved simultaneously:
with
,$,y,., =1.oo;
equations
(12) Condenser
for the reboiler,
4 Y, = K,X,
-
(13)
(18)
+ l)+An,.A,,
dX, --=!!!(Y,-X,)+R,-X,.An,, dt A,
with
i$Yhi =
= I’, -D(R
dt
1.Oo;
(19)
(14)
for the condenser,
Y, = K,X,
+Qo-fd$
(15)
with
0 = h,x,, -
o= Constraint
l.OO-
Y,,,
(20)
3
1 i =
1,4,
(21)
cyo,, i= I
(22)
on the enthalpy
(23)
If the temperature and the composition of the liquid phase are known, then the molar enthalpy of this liquid phase is determined. At each time, the molar enthalpy calculated according to the value of the composition (given by the integration of the mass-balance equations) and the temperature (given by the liquid-vapor equilibrium constraint equations) must be equal to the molar enthalpy calculated by the integration of the energy-balance equation: HI,, = Ifi (Xl> T,, pj). Total enthalpies heat of reaction (6) and (9).
and
Plate j, 1
dX,_ 1
(17)
dt -A,
are used in this formulation; thus, no terms are required in equations (3)
NUMERICAL
v,+,(Y,+,
- V,(Y,-X,)
(25)
-Xj)+L,ml
(26)
1
(27)
+R,-X;An,,
-L,-,H,,-,
SOLUTION
The method proposed by Distefano [2] and also adopted by Boston et al. [4] used the following approximations in order to reduce the system to one in which the differential equations can be integrated separately from the algebraic equations. The calculation of the liquid-vapor equilibrium is performed outside of the integration. That is, the equilibrium temperature and the composition of the vapor phase are calculated for a given composition of the liquid phase and pressure (bubble-point calculation). With this temperature, the molar enthalpy is calculated, using the constraint equation on the enthalpy. The derivative of the enthalpy with respect to time is next estimated using a low-order approximation (e.g. finite difference). Then the enthalpybalance equation is solved separately as an algebraic equation. With this approach, only the mass-balance equations need to be integrated. In batch rectification with chemical reactions, this kind of approximation is inappropriate: the variation of the enthalpies will be too rapid to be estimated with acceptable precision. Therefore, the following
(24)
0 = A, - A (X,, J’,, T,,, r,,);
- YH,,
(28)
-L,H,;,+Q,-4,~
0 = y, + ,,,( 1 - e,) + e,k,,x,: - Y,~, i = 1, ncr (29)
o=
l.OO-
c yjj, r=l
(30)
0 = H,j - H,(X,, P,, 7’J
(31)
0 = A, - A (X,, P,, 7;, 7,);
(32)
and
Reboiler dA, =L,dt
V,+F+An,.A,,
dX, 1 = 7 MX” - X,) - V,(Y, - X,), dt ,, +F(Xr--XJ
1
+R,--Xh.Anh,
(33)
(34)
(35)
P. E. CUILLEand G. V.
392
dH,, 2 dt
1 =A,
the initial column conditions and a device for implicitly generating initial derivative values of the algebraic variables.
L,,H,,” - VbHv,b + FH, dA,
+Qdf,,t,~
>
O=k,,x,,--y,,,
o= l.OO-
>
Calculation
i = l,n,,
(37)
i y,, r=l
(38)
and 0 = H,.b - H, (X,, Pi,, r,).
(39)
It should be noted that since both the total mole and the species mole balances are used, the condition requiring the sum of the mole fractions on a plate to be equal to one is redundant. In fact, it is easy to show that the sum of each of the n, time derivatives [equations (19) (26) and (34)] must always be equal to zero. Thus, the sum of the mole fractions will always be equal to its initial value of one. Numerical
method
This set of differential and algebraic equations is solved using the routine “LSODI” from the package “ODEPACK” created at the Lawrence Livermore Laboratory [7]. This program uses either Gear’s method (also known as the backward differentiation formula or BDF) or the implicit Adam’s method. The BDF method is numerically the most efficient method for stiff systems: on a given simulation example the BDF is three times faster than Adam’s method. This program solves systems of the following form:
A.(t, Y)* % = G(t, Y), where A\ is a square matrix, Y and G are vectors and t the independent variable (e.g. time). If the system contains at least one algebraic equation, the matrix A will be singular. CALCULATION
OF THE INITIAL
REKLAITIS
CONDITIONS
Before starting the integration of the balance equations, it is necessary to calculate a consistent set of initial values for all the variables of the system (for both algebraic and differential equations). Moreover, because of the algebraic equations, the integration program is not able to calculate all the initial time derivatives because the A matrix is singular. Therefore, it is necessary to devise a procedure for calculating those derivatives, before being able to start the integration. In previous work, this problem is treated by setting all initial derivatives of the algebraic variables to zero, choosing a small initial integration step size, and bypassing the integrator error tolerance checks for the first integration step [6]. This ad hoc procedure is not reliable in general and is especially inappropriate in the presence of reactions. In this section we discuss both a methodology for calculating
of the initial vaIues
We consider, as a starting point, the column operating at total reflux in a steady state and without any chemical reactions. Consequently the problem of determining the initial values is equivalent to calculating this steady state, with the same assumptions as those used in the dynamic model (negligible vaporphase holdup, constant holdup on each plate). When the column is operating at total reflux, the holdup and the composition in the reboiler are constants with time. The reboiler mass-balance equations, therefore, reduce to L, = Vb
(41)
x, = Y,.
(42)
and
If we consider that the total initial charge is in the reboiler, then Y, corresponds to the bubble-point composition of the initial charge. By virtue of the above equality, the liquid composition X, is known and thus from a bubble-point calculation Y, will be obtainable. Again, since the plate holdup and composition are constant with time, the plate mass-balance equations for plate j reduce to L,=
v,+i
(43)
and X,=Y,+,.
(44)
Thus, with Y, known, X,_ , is determined and from a bubble-point calculation Y,_ , is determined. Consequently, by proceeding from stage-to-stage alternating bubble-point calculations with use of the trivial equality X, = Y, + ,, the initial estimate of the column composition profile can be found. Using this composition profile and the bubble temperatures, the liquid densities on each plate can be calculated and, since the plate holdup volume is specified, the molar holdup of each component can be determined. Overall mass balances will then establish a new composition and a new holdup for the reboiler. If necessary we repeat the calculation of the composition profile until convergence is attained. The convergence conditions are based on the variations in the reboiler holdup and compositions between two successive iterations. Once the column composition and holdup calculations have converged, the flow rates between the plates can be calculated by using the steady-state form of the energy-balance equations (28) together with equation (44). Because of conditions [44] the plate energy balances are uncoupled and this calculation simply amounts to a straightforward sequential solution of the stage energy-balance
393
Simulation of a staged batch distillation unit equations beginning with either the reboiler or the condenser. The complete calculation essentially amounts to a successive substitution scheme in which the reboiler holdup and composition are the iteration variables. If the initial charge of the reboiler is large compared with the total holdup of all the plates, the holdup and composition in the reboiler at total reflux are very close to the initial charge, and the calculation converges very quickly. But if the initial charge is comparable with the total holdup of the plates, the calculation may be nonconvergent: stable oscillations around the solution may occur. This type of numerical behavior is well-known in equation-solving applications and the standard solution is the use of relaxation factors. This, however, requires the empirical tuning of relaxation factors which typically have little physical significance. In the present instance, a physically interpretable relaxation procedure can be formulated. Essentially the numerical oscillations arise because the algorithm instantaneously fills up the plates and thus, in cases with small initial charge, significantly alters the reboiler holdup and composition. Physically, of course, this would not happen: the trays would fill gradually. A reasonable and very stable relaxation procedure is to only fill a fraction of the total tray holdup with each update of the composition. The fraction can be selected by comparing the total specified tray molar holdup (using the density of the initial charge) with the reboiler charge. If ZZ:A,/A,is small (say, < 0.10) no relaxation is required. Otherwise, a generally safe relaxation factor, in our experience, is about l/2, corresponding to filling half of the yet unfilled tray holdup with each composition iteration. In summary, the initial reflux calculation proceeds as follows: Step 0. Set A, = initial charge X, = initial charge composition A to)= 0
I
In our experience this has proven to be a very reliable and quick method. By way of example, Fig. 2 shows the variation of the value of the reboiler holdup (in moles) vs the number of iterations, with and without the underrelaxation (respectively A = 0.5 and 1 = 1). Calculation
of the initial derivatives
The system
to be solved
The elements
of the A matrix
are
if i #j, if the ith equation
is differential,
a,j = 0,
if the jth equation
is algebraic.
and
In order to numerically solve the system of equations (45) beginning at time t =O, it is necessary to use values of dY/dt at t = 0. The initial derivative values can, of course, be readily obtained for the yi corresponding to the differential equations by merely evaluating the associated g,(Y). For the algebraic equations, the evaluation is problematic. One approach is to observe that for the algebraic variables y,, since g,(Y) = 0, dg, = z and from the chain function of t,
dt = 0
rule, since g, is not an explicit
Thus,
$50 Y, Set
Step 2. For plates j = n through 0: (a) Compute the bubble compositon Y, and temperature in plate j. (b)SetX,_,=Y,andj=j-1 andgo to (a). Step 3. (a) Compute the molar holdup A:‘). (b) Check for convergence. If converged, go to Step 5. (c) Set A(“)=A”~‘+(A!“)-_A(“-‘) / I J J 1.1 calStep 4. Perform the overall mass-balance culations to update A, and X, and go to Step 1. using
(45)
a,=O, a, = 1,
'dt
Step 1. Compute the bubble composition and temperature of the reboiler. x, = Y,.
form:
A.%=G(Y).
.
the stage flows Step 5. Calculate stage energy balances.
has the following
the
where
’
J, is the ith row of the Jacobian
matrix
of
G W. 500
=
450
E 4 12 0 400 .c &I .C :: ($
350 Relax * 1 3oc I
III1 01234
1 5
II
II 6
7
6
1IIlI 9
10
11
12
13
14
Iterations Fig. 2. Effect of relaxation factor on the calculation initial conditions.
of the
P. E. CUILLE and G. V. REKLAITIS
394 The above equations trivial equations,
can be combined
dyp
with the
dy,
dt=dt> for the differential
variables
to obtain
the system
G*=J*dy. dt
where J* is a matrix whose ith row consists of the coordinate vector eT if the ith equation is differential and of J,, if it is algebraic; and where, g: = dyi’)/dt if the ith equation is differential and g: = 0, if thejth equation is algebraic. This relation is a linear algebraic system, in the vector of unknowns dY/dt. If the J* matrix is nonsingular, the solution of this system will give all the initial derivatives. In the case of the system of equations (10)(29), the matrix J* is singular (the system has an infinity of solutions, including dY/dt). One method for solving this problem is to estimate a subset of derivatives in order to reduce the dimension of the system and to obtain a resulting nonsingular system. For instance, it is sometimes possible to assume that some derivatives are equal to zero at the initial time. However, this method is ad hoc and can not be guaranteed to work in all cases. Moreover, this kind of approximation is not possible with chemical reaction since with chemical reactions the initial derivatives will, in general, not be zero. Therefore, the following approach is preferable. The mathematical cause of this numerical problem is the discontinuity of the time derivatives of the solution functions at the starting point. This discontinuity has no physical reality since the start of production corresponds to the opening of a valve, which physically requires a small but nonzero duration. We can simulate the valve opening, by using a continuous function for the corresponding flow rate which increases the flow from zero to a finite positive value, in a small but nonzero time interval. The time derivatives of this function will be equal to zero at the beginning and at the end of this interval. For instance, if the column is to be simulated with constant production rate D = Do, then over the arbitrary time interval from, say, t = - 1 s to t = 0 the variable D in equation (18) can be replaced with the cubic polynomial function: D = D,( -2t3
- 3t2 + 1);
at
t=
-1,
D=O
and
g=O dt
and
t =O,
D =D,
and
dD -=O. dt
At the beginning of this time interval, i.e. at t = - 1, the column is operating at total reflux: all the time derivatives are equal to zero and are continuous functions.
Over the time interval - 1 to 0, the model equations are integrated with the cubic forcing function for D inserted into the model. At time t = 0 appropriate values of all derivatives of algebraic variables will have been calculated through the normal execution of the integrator and, thus, the integration can continue without interruption. In this way both the integration at t = - 1 s (total reflux condition) and at t = 0 (beginning of batch production) can be initiated without any numerical difficulties and without explicit estimates of the derivatives of algebraic variables. It should further be noted that this device can also be employed if other column specifications are of interest. For instance, if the column heating rate is specified as fixed and the distillate rate is left variable, then the reflux ratio variable can be subjected to the cubic forcing function. Specifically, at t = - 1, l/R = 0 since the column is at total reflux and at t = 0, l/R will take on the specified value. The transition function for l/R will thus correspond to the opening of a valve activated by a proportional controller. The only minor error introduced using this transition device is the fact that a distillate quantity equal to the integral of the distillate production curve over the time interval - 1 to 0 will have been produced prior to time 0. However, since the transition time interval can always be chosen to be small compared to the total batch distillation time, this error can be made as small as desired. (In our examples, we used 1 s for a 9000 s operating time.) EXAMPLES
OF SIMULATION
A general program simulating a batch rectification column has been developed, using the model and the numerical methods presented above. The program runs on a Control Data “Cyber 205” supercomputer and is written in FORTRAN 77 extended for vector calculations. This program can be easily adapted to any physical properties package, using small interface subroutines. The regular version uses the physical properties package of the simulator FLOWTRAN [8] through an interface which converts American units to SI units and adds heats of formation to the predicted enthalpies when reacting systems are indicated. Several examples of rectification with and without chemical reactions have been simulated, and will be reported next. In the first two examples, the options of the FLOWTRAN physical properties package used are the following: Antoine vapor pressure; Redlich-Kwong vapor fugacity; corrected liquid fugacity; ideal solution activity coefficients. Example
I
We have simulated successfully the rectifications performed by Domenech and Enjalbert [9] involving
Simulation of a staged batch distillation unit
Production
(mol)
Fig. 3. Simulation Example I. Experimental points and simulated curve for R = 3 and x,,, = 0.62 for C,H,,.
9000
Time
(s)
Fig. 6. Simulation Example time and plates for
1. Mole fraction
of C,H,,
vs
R = 2 and x,,, = 0.55.
the separation of cyclohexane and toluene. The experimental equipment used was a perforated-plate column, with 4-10 plates with a 60 1. reboiler heated with a heat transfer rate of 3 kW; this column operates under atmospheric pressure with a reflux ratio between 3 and 10 (see Ref. [9] for further details of the operating conditions). Figures 3-5 show the results of the simulations (curves) and the experimental data with reflux ratios of 3, 4 and 6 for different initial composition (respectively: molar fraction in cyclohexane of 0.62,0.55 and 0.42). Figures 68 show the variations of the
1
9000
o-
09- tb
Time
uou
080.7
6 06-
(s)
Fig. 7. Simulation Example time and plates for
1. Mole fraction
of C,H,,
vs
R = 4 and x,,, = 0.55.
-
E 058 G
04-
x” 0302Ol00
I 0
100
50
Production
150
200
(mol)
Fig. 4. Simulation Example 1. Experimental points and simulated curve for R = 4 and x,,, = 0.55 for C,H,,.
Time
(sl
Fig. 8. Simulation Example time and plates for
ooi-----k
Production
I 100
I 150
(mol)
Fig. 5. Simulation Example 1. Experimental points and simulated curve for R ==6 and x,,~ = 0.42 for C,H,,.
1. Mole fraction
of C,H,,
vs
R = 6 and x,,, = 0.55.
composition profile of cyclohexane in the column during the simulation, on a given example with reflux ratios of 2,4 and 6 (with the same initial cyclohexane molar fraction: 0.55). The flat portion of the surfaces drawn in Figs 6-8 corresponds to a zone in which the column performs a good separation. It is evident that the size of this flat zone increases with the reflux ratio, but this increase is achieved at the cost of higher utilities requirements.
P. E. CUILLE and G. V. REKLAITIS
396
Example 2 We also tested simulation cases with chemical reactions: for example the esterification of acetic acid with I-propanol which also involves the undesirable hydrolysis reaction: CH,-COOH
+ CH,--CH+ZH,OH
+ CH,-CQO--CH2-CH,-CH,
+ H,O
(46)
and CH,-COO-CH,-CH,--CH, + CH,-COOH
+ H,O
+ CH,-CH,-CH,OH.
(47)
We calculate the initial state at total reflux with a mixture of acetic acid and water (therefore without chemical reaction). Then we start the simulation in a fed batch mode: we feed the reboiler with a mixture of propanol and water and simulate the rectification with the production of propyl acetate. Figure 9 shows the molar fraction of ester vs the plate numbers and time. At t = 0, there is no ester in the column. In the reboiler [Plate 6) the reactions start immediately, but at the top of the column (condenser = Plate 0) the start of reaction is delayed. In the top plates, where the molar fraction of water is important, we can see that the fraction of ester vs time reaches a maximum, and then decrease, because of the hydrolysis reaction. Figure 10 shows the variation of the molar fraction of water. At t = 0, the column contains only acid and water which is the most volatile compound and the fraction of water is maximal at the top of the column. Just after the start of the rectification, because of the propagation of alcohol in the column, the fraction of water on the top plates falls down very quickly. Later, the fraction of water becomes maximum on Plate 2 because the alcohol is now the most volatile compound.
gooo sw Fig. 10. Simulation
Example 2. Mole fraction time and plates.
The initial parameters
of water
for this simulation
vs
are:
stages: reboiler, 5 plates, condenser; holdup of the condenser and of each plates = 0.25 dm3; Murphree efficiency of the plates = 70%; reflux ratio = 3; top pressure = 101,300 Pa, no pressure drop between plates; initial charge = 200 mol (acid 90%, water 10% mole fractions); composition of the feed = alcohol 90%, water 10% (mole fractions); feed rate = 0.02 mol SK’; esterification rate constant = 2 x hydrolysis rate constant (independent of the temperature). Example 3
p 5
Fig. 9. Simulation
Example 2. Mole fraction and plates.
of ester vs time
As a further illustration, we present a case involving the linkage of the simulation model to the more extensive thermodynamic properties package described by LeLann et al. [lo, 111. The simulation was executed on a Bull SPS9 computer operating under UNIX. The example involves the separation of the ethanol-water azeotrope from an initial mixture containing 10% (mol) ethanol in water. For this application, the liquid-phase behavior is modeled using the UNIQUAC method and the vapor-phase behavior using the Peng-Robinson equation of state. The results (not shown here) confirm that for operation of the column at 1 atm, one can, of course, obtain reasonably accurate results using the ideal gas equation of state. However, near the azeotropic condition, Peng-Robinson is clearly more accurate. Figure 11 shows the liquid composition profile for ethanol at the initiation of rectification and after 30 min of production in a 20-plate column. As might be expected, the initial flat profile tending towards the
391
Simulation of a staged batch distillation unit 1
.o
r
09 0.8 07 06 '05 0.4 0.3 02 0.1 0.0 0
1 2 3 4 5 6
7 8 9101112131415161718192021 Plates
Fig. 11. Simulation Example 3. Mole fraction of ethanol vs plates.
azeotropic composition becomes steeper over time as the alcohol-rich distillate is removed. USE OF THE “CYBER
205” SUPERCOMPUTER
The main advantage of this type of supercomputer, compared with the classical “Von Neuman” machine, is the capability to perform parallel calculations when executing operations on vectors. On the Cyber 205, two ways of using this capability are available:
I. The automatic
“vectorization” of the do-loops: if a do-loop satisfies certain vectorization constraints, the Fortran compiler can translate it into a machine vector instruction. If those vectorization constraints are satisfied, programs written in standard Fortran 77 (or Fortran 66) can use the capabilities of the vector processor. 2. The explicit vector programming: using some extensions of the Fortran language (specific to the Cyber 205) it is possible to program directly in a vectorized form. The integration package used in the simulation was written in standard Fortran 77, however, it is not well-adapted to the vector capabilities of the machine. Since only a few of the do-loops are vectorizable, the use of the automatic vectorization of the loops on a given example of simulation (the cyclohexane-toluene separation) saves only 3% of the CPU time compared with a simulation without this vectorization. In order to enhance the use of the vector capabilities, we rewrote the solver used by the integrator at the corrector step. The original one used a Gaussain elimination method (GEM) and made very poor use of the vector processor. The new one uses a Jordan method which is numerically less efficient on a classical computer, but is better adapted to be programmed in a vectorial form. For instance, on a separate example involving a system of 50 algebraic linear equations, this routine is more than four times faster than that based on a GEM. Using this routine in the simulation, we can save more than 40% of the CPU time on the cyclohexane-toluene example.
The actual version of the simulation program contains 2081 lines of Fortran code (without the integrator and the FLOWTRAN physical properties package). The CPU time needed for a simulation depends on the characteristics of the simulated problem. For the same column (5 plates) the program need only 3 s of CPU for simulating 9000 s of operation for the conditions of Example 1 (two compounds, no chemical reaction), and 73 s of CPU for the conditions of Example 2 (four compounds, two chemical reactions) also for 9000 s of operation. A better use of the vector processor is limited by a major bottleneck: the calculation of the physical properties of the gas and liquid mixtures used in the rectification. There is no physical properties package written in a vectorized form available today. If the use of computers with vector capabilities become more widespread, it will be necessary to develop new packages able to calculate not only one physical property at a time, but also a series of values of a given property for a set of different mixtures. For example: to calculate simultaneously all the liquid enthalpies for all the liquid mixtures present in the rectification column at a given time. Acknowledgemenrs-This work was carried out with the financial support of Rhone Poulenc, S.A., which is gratefully acknowledged. The work was completed at Purdue University while the first author was on leave as a visiting scholar and is part of the research presented in his Docteur Thesis.
NOMENCLATURE
A D e F H, H,
= = = = = =
K= k = L = M = n = n, = neq = P = $ 1
Material holdup (mol) Overhead product flow rate (mol s-‘) Murphree efficiency Feed rate Molar enthalpy of the liquid phase Molar enthalpy of the vapor phase Vector of k-values Element of the vector K Liquid flow rate (mol s-‘) Molar mass Number of plates Number of components Number of equations Pressure (Pa) ;;Fu;frit;tt transfer (W)
R = Vector of chemical reaction rates s = Iteration counter T = Thermodynamic temperature (K) t = Time(s) V = Vapor flow rate (mol s-l) X = Vector of mole fractions (liquid) x = Element of the vector X Y = Vector of mole fractions (vapor) y = Element of the vector Y An = ZR, r = Volume (m’) p = Density (kg mm3) 1 = Relaxation factor Subscripts 0 = Condenser 1 = Plate 1 (top) j = Plate j
P. E. CUILLE and G. V. REKLAITIS
398 n b i ini f
= = = = =
Plate n (bottom) Reboiler Component i Initial Feed
REFERENCES 1. E. L. Meadow, Chem. Engng Prog. Symp. Ser. 46(59), 48 (1963). 2. G. P. Distefano, AIChE JI 14, 190 (1968). 3. C. W. Gear. Communs Ass. comout. Mach. 14. 176 (1971). 4. J. F. Boston, H. I. Britt, S. Jirapongphan and V. B. Shah, An Advanced System for the Simulation of Batch
Distillation Operation, Vol. 11, pp. 203-237. Foundation of Computer Aided Chemical Process Design, Engineering Foundation, New York (1980). G. D. Byrne, A. C. Hindmarsh and K. C. Jackson, _J. Comput. ihem. Engng 1, 133 (1977). 6. S. E. Gallun and C. D. Holland, Comput. them. Engng 6, 231 (1982). Ass. comput. Mach. Signum News1 15, 7. A. C. Hindmarsh, 10 (1980). 8. J. D. Seader, W. D. Seider and A. C. Pauls, FLOWTRAN Simulation, an Introduction, 2nd edn. CACHE, Ann Arbor, Mich. (1977). 9. S. Domenech and M. Enjalbert, Chem. Engng Sci. 29, 1519 (1974). Docteur Ingenieur Thesis, INP 10. J. M. LeLann, Toulouse, France (1983). 11. J. M. LeLann, X. Joulia and B. Koehret, Entropic 123, 29 (1985).
Copynght 0
0098-l 354/86 $3.00 + 0.00 1986 Pergamon Journals Ltd
DYNAMIC SIMULATION OF MULTICOMPONENT BATCH RECTIFICATION WITH CHEMICAL REACTIONS P. E. CUILLE Institut du Genie
Chimique,
Institut
National
Polytechnique
de Toulouse,
Toulouse,
France
G. V. REKLAITIS School
of Chemical
Engineering,
Purdue
University.
West Lafayette,
IN 47907,
U.S.A.
(Received 6 February 1985; revision received 29 October 1985; received,for publication 20 December 1985) Abstract-A model is formulated and solution strategies are investigated for the simulation of a staged batch distillation unit with chemical reactions in the liquid phase. The implementation of the model involves three computation phases: calculation of column profiles at total reflux, calculation of initial derivatives of the algebraic variables and integration of the dynamic model itself. The latter calculation is performed by solving the system of algebraic and differential equations simultaneously using an implicit integration method. The possibilities of vectorizing the model are considered and results are presented on the use of parallel computation. Examples with and without reactions on the plates are reported. Scope-For centuries, batch processes were used in a very wide field of application, from food processing to coal and petroleum engineering. Beginning at about World War II. as it became necessary to cheaply and quickly produce large quantities of products, continuous processes became the standard. More recently, because of the development of new processes for the production of small amounts of products with high added value, batch processes returned to fashion. Publications about the mathematical modeling of batch distillation operations without chemical reactions are found in the literature, as early as the sixties. One of the first rigorous mathematical models was proposed by Meadow [l]. The solution of such a model, using various numerical methods, has been presented by Distefano [2]. During the same period Gear [3] published his first work on a numerical method for the integration of ordinary differential equations. More recently, because of the increasing interest in batch processes, some new works were published: Boston et al. [4] proposed the program “BATCHFRAC” as a general simulation package for batch rectifications. This program uses the model proposed by Distefano, with some numerical approximations, and employs the integration routine “EPISODE” [5], based on the numerical method proposed by Gear. The key approximation employed in Boston’s work consists of the calculation of the time derivative of the enthalpies via a low-order difference approximation. Gallun and Holland [6] proposed a more rigorous approach, using contemporary improvements in the numerical solution of sets of coupled algebraic and differential equations, resulting from Gear’s works. Gallun’s approach is thus based on the direct solution of the balance equations (differential) and of the constraint equations (algebraic) describing the system, He also uses a hydrodynamic model for the calculation of the pressure drop between each plate and of the volumetric holdup on each plate. The plates are assumed to be ideal stages for the vapor-liquid equilibrium calculation. Reactions on the plates are not included in this model. In the present case we do consider reactions and nonideal stages, and also examine the possibilities for vectorization. However, neither this work nor the previously cited references consider the dynamics of the column heat-up or start-up phases which result in the attainment of the total reflux condition. Conclusions and Significance-The simultaneous solution of the model algebraic and differential equations by means of the Implicit integrator LSODI proves to be a highly effective and reliable approach for batch rectification with reaction. Calculation of the initial, total reflux profile can be time consuming if the initial reboiler charge is not large compared to the holdup on the plates. However, a simple physically meaningful relaxation approach provides satisfactory convergence in all test cases examined in this work. The problem of estimating the initial derivatives of the algebraic variables could be very reliably resolved by introducing a staightforward polynomial transition function. In general, the computation time required to simulate a column of a given number of stages over a given simulation time interval increases significantly but manageably with number of components and number of reactions. The parallel computation trials showed that the principal benefits of this mode of computation are achieved if the linear equation solver used in the corrector step is recoded in vectorized form. With such modifications, some 40% can be saved in execution time. Use of automatic vectorization yielded only minor execution time savings (3%). Further significant gains due to vectorization can only be expected if properties evaluations can be carried out in vector-compatible form.
THEORETICAL
MODEL
In this section we consider the mass- and energybalance equations for the batch rectification column
consists of the total mass balance and mass balance on each component. The energy-balance differential equation merely consists of the applica-
shown
tion
in Fig.
1. The
set of mass-balance
equations
differential 389
of
the
first
law
of
thermodynamics.
The
390
P. E. CUILLE and G.
V.
REKLAITIS
and dA,.H,, L dt
f-4
(3)
Plate j, 1
A0
*
I
dJ&= dt
w
0
‘0
l_
= V,.H,,,-(L,+D).H,,~+Q,;
(4)
V,+l'Y,+l+Lj-,'X,~,-V,Y, -
L;X,+R;A,
(5)
and w=
V/+I
dt
i
0,
. ffv., + I +Ljm,‘ffl,,m,
- v,.Hv,,--,.4.,$-Q,;
(6)
Reboiler dA, = L,-dt dA,.X, ~ dt
T/,+An,.A,,
=L;X,-
(7)
V,,.Y,+R,.A,
(8)
and dA,. H1.b ~ = L;H,,” dt Fig. 1. Column schematic.
equations following
presented hereafter assumptions:
are valid
under
the
1. The vapor-phase holdup is assumed to be negligible compared to the liquid-phase holdup on each plate. 2. Chemical reactions in the vapor phase will be neglected. 3. The initial state of the column is the steadystate total reflux condition with no reactions. 4. The liquid volumetric holdups on the plates will be assumed to be constant. Thus, our model is directed at simulating the dynamics of the main production period during which the hydrodynamic conditions are not widely varying. 5. The pressure drops and the plate efficiencies are constant during the operation. The following
differential
balance
equations
result.
Condenser dA, -=V,-L,--D+An,.A, dt = V,An,.A,--D(R dA,.X, -= dt
V,.Y,-L,.X,-D.X,+R,.A,
(1) (2)
(9)
+ Q,,.
An analysis of the degree of freedom of this system readily shows that with these differential equations alone the problem is underspecified. In order to obtain a fully determined system, the variables appearing in these balance equations must also satisfy the following equations. Constraint
on the volume
The volume of the liquid phase in the condenser and on each plate is assumed to be constant (we do not use a hydrodynamic model). Consequently the corresponding molar holdups are functions only of the temperature, the pressure and the composition: A, = r, P,IM,
and
P/ = P CX,, T,, f’,).
(10)
This holdup must be equal to that calculated by the integration of the total mass-balance equations. It should be noted in passing that the inclusion of tray hydrodynamic models such as the Francis-Weir equation merely results in additional equations which can be appended at this junction [e.g. 61. The numerical solution approach discussed in the next section will not be materially altered. Constraint
+ l),
- I”,.&,
due to vapor-liquid
equilibria
For given values of the temperature, and the composition: for each plate j, Y,=Y,+,.(l
-e,)+e;K;X,
the pressure
(11)
391
Simulation of a staged batch distillation unit system of linked differential and algebraic has to be solved simultaneously:
with
,$,y,., =1.oo;
equations
(12) Condenser
for the reboiler,
4 Y, = K,X,
-
(13)
(18)
+ l)+An,.A,,
dX, --=!!!(Y,-X,)+R,-X,.An,, dt A,
with
i$Yhi =
= I’, -D(R
dt
1.Oo;
(19)
(14)
for the condenser,
Y, = K,X,
+Qo-fd$
(15)
with
0 = h,x,, -
o= Constraint
l.OO-
Y,,,
(20)
3
1 i =
1,4,
(21)
cyo,, i= I
(22)
on the enthalpy
(23)
If the temperature and the composition of the liquid phase are known, then the molar enthalpy of this liquid phase is determined. At each time, the molar enthalpy calculated according to the value of the composition (given by the integration of the mass-balance equations) and the temperature (given by the liquid-vapor equilibrium constraint equations) must be equal to the molar enthalpy calculated by the integration of the energy-balance equation: HI,, = Ifi (Xl> T,, pj). Total enthalpies heat of reaction (6) and (9).
and
Plate j, 1
dX,_ 1
(17)
dt -A,
are used in this formulation; thus, no terms are required in equations (3)
NUMERICAL
v,+,(Y,+,
- V,(Y,-X,)
(25)
-Xj)+L,ml
(26)
1
(27)
+R,-X;An,,
-L,-,H,,-,
SOLUTION
The method proposed by Distefano [2] and also adopted by Boston et al. [4] used the following approximations in order to reduce the system to one in which the differential equations can be integrated separately from the algebraic equations. The calculation of the liquid-vapor equilibrium is performed outside of the integration. That is, the equilibrium temperature and the composition of the vapor phase are calculated for a given composition of the liquid phase and pressure (bubble-point calculation). With this temperature, the molar enthalpy is calculated, using the constraint equation on the enthalpy. The derivative of the enthalpy with respect to time is next estimated using a low-order approximation (e.g. finite difference). Then the enthalpybalance equation is solved separately as an algebraic equation. With this approach, only the mass-balance equations need to be integrated. In batch rectification with chemical reactions, this kind of approximation is inappropriate: the variation of the enthalpies will be too rapid to be estimated with acceptable precision. Therefore, the following
(24)
0 = A, - A (X,, J’,, T,,, r,,);
- YH,,
(28)
-L,H,;,+Q,-4,~
0 = y, + ,,,( 1 - e,) + e,k,,x,: - Y,~, i = 1, ncr (29)
o=
l.OO-
c yjj, r=l
(30)
0 = H,j - H,(X,, P,, 7’J
(31)
0 = A, - A (X,, P,, 7;, 7,);
(32)
and
Reboiler dA, =L,dt
V,+F+An,.A,,
dX, 1 = 7 MX” - X,) - V,(Y, - X,), dt ,, +F(Xr--XJ
1
+R,--Xh.Anh,
(33)
(34)
(35)
P. E. CUILLEand G. V.
392
dH,, 2 dt
1 =A,
the initial column conditions and a device for implicitly generating initial derivative values of the algebraic variables.
L,,H,,” - VbHv,b + FH, dA,
+Qdf,,t,~
>
O=k,,x,,--y,,,
o= l.OO-
>
Calculation
i = l,n,,
(37)
i y,, r=l
(38)
and 0 = H,.b - H, (X,, Pi,, r,).
(39)
It should be noted that since both the total mole and the species mole balances are used, the condition requiring the sum of the mole fractions on a plate to be equal to one is redundant. In fact, it is easy to show that the sum of each of the n, time derivatives [equations (19) (26) and (34)] must always be equal to zero. Thus, the sum of the mole fractions will always be equal to its initial value of one. Numerical
method
This set of differential and algebraic equations is solved using the routine “LSODI” from the package “ODEPACK” created at the Lawrence Livermore Laboratory [7]. This program uses either Gear’s method (also known as the backward differentiation formula or BDF) or the implicit Adam’s method. The BDF method is numerically the most efficient method for stiff systems: on a given simulation example the BDF is three times faster than Adam’s method. This program solves systems of the following form:
A.(t, Y)* % = G(t, Y), where A\ is a square matrix, Y and G are vectors and t the independent variable (e.g. time). If the system contains at least one algebraic equation, the matrix A will be singular. CALCULATION
OF THE INITIAL
REKLAITIS
CONDITIONS
Before starting the integration of the balance equations, it is necessary to calculate a consistent set of initial values for all the variables of the system (for both algebraic and differential equations). Moreover, because of the algebraic equations, the integration program is not able to calculate all the initial time derivatives because the A matrix is singular. Therefore, it is necessary to devise a procedure for calculating those derivatives, before being able to start the integration. In previous work, this problem is treated by setting all initial derivatives of the algebraic variables to zero, choosing a small initial integration step size, and bypassing the integrator error tolerance checks for the first integration step [6]. This ad hoc procedure is not reliable in general and is especially inappropriate in the presence of reactions. In this section we discuss both a methodology for calculating
of the initial vaIues
We consider, as a starting point, the column operating at total reflux in a steady state and without any chemical reactions. Consequently the problem of determining the initial values is equivalent to calculating this steady state, with the same assumptions as those used in the dynamic model (negligible vaporphase holdup, constant holdup on each plate). When the column is operating at total reflux, the holdup and the composition in the reboiler are constants with time. The reboiler mass-balance equations, therefore, reduce to L, = Vb
(41)
x, = Y,.
(42)
and
If we consider that the total initial charge is in the reboiler, then Y, corresponds to the bubble-point composition of the initial charge. By virtue of the above equality, the liquid composition X, is known and thus from a bubble-point calculation Y, will be obtainable. Again, since the plate holdup and composition are constant with time, the plate mass-balance equations for plate j reduce to L,=
v,+i
(43)
and X,=Y,+,.
(44)
Thus, with Y, known, X,_ , is determined and from a bubble-point calculation Y,_ , is determined. Consequently, by proceeding from stage-to-stage alternating bubble-point calculations with use of the trivial equality X, = Y, + ,, the initial estimate of the column composition profile can be found. Using this composition profile and the bubble temperatures, the liquid densities on each plate can be calculated and, since the plate holdup volume is specified, the molar holdup of each component can be determined. Overall mass balances will then establish a new composition and a new holdup for the reboiler. If necessary we repeat the calculation of the composition profile until convergence is attained. The convergence conditions are based on the variations in the reboiler holdup and compositions between two successive iterations. Once the column composition and holdup calculations have converged, the flow rates between the plates can be calculated by using the steady-state form of the energy-balance equations (28) together with equation (44). Because of conditions [44] the plate energy balances are uncoupled and this calculation simply amounts to a straightforward sequential solution of the stage energy-balance
393
Simulation of a staged batch distillation unit equations beginning with either the reboiler or the condenser. The complete calculation essentially amounts to a successive substitution scheme in which the reboiler holdup and composition are the iteration variables. If the initial charge of the reboiler is large compared with the total holdup of all the plates, the holdup and composition in the reboiler at total reflux are very close to the initial charge, and the calculation converges very quickly. But if the initial charge is comparable with the total holdup of the plates, the calculation may be nonconvergent: stable oscillations around the solution may occur. This type of numerical behavior is well-known in equation-solving applications and the standard solution is the use of relaxation factors. This, however, requires the empirical tuning of relaxation factors which typically have little physical significance. In the present instance, a physically interpretable relaxation procedure can be formulated. Essentially the numerical oscillations arise because the algorithm instantaneously fills up the plates and thus, in cases with small initial charge, significantly alters the reboiler holdup and composition. Physically, of course, this would not happen: the trays would fill gradually. A reasonable and very stable relaxation procedure is to only fill a fraction of the total tray holdup with each update of the composition. The fraction can be selected by comparing the total specified tray molar holdup (using the density of the initial charge) with the reboiler charge. If ZZ:A,/A,is small (say, < 0.10) no relaxation is required. Otherwise, a generally safe relaxation factor, in our experience, is about l/2, corresponding to filling half of the yet unfilled tray holdup with each composition iteration. In summary, the initial reflux calculation proceeds as follows: Step 0. Set A, = initial charge X, = initial charge composition A to)= 0
I
In our experience this has proven to be a very reliable and quick method. By way of example, Fig. 2 shows the variation of the value of the reboiler holdup (in moles) vs the number of iterations, with and without the underrelaxation (respectively A = 0.5 and 1 = 1). Calculation
of the initial derivatives
The system
to be solved
The elements
of the A matrix
are
if i #j, if the ith equation
is differential,
a,j = 0,
if the jth equation
is algebraic.
and
In order to numerically solve the system of equations (45) beginning at time t =O, it is necessary to use values of dY/dt at t = 0. The initial derivative values can, of course, be readily obtained for the yi corresponding to the differential equations by merely evaluating the associated g,(Y). For the algebraic equations, the evaluation is problematic. One approach is to observe that for the algebraic variables y,, since g,(Y) = 0, dg, = z and from the chain function of t,
dt = 0
rule, since g, is not an explicit
Thus,
$50 Y, Set
Step 2. For plates j = n through 0: (a) Compute the bubble compositon Y, and temperature in plate j. (b)SetX,_,=Y,andj=j-1 andgo to (a). Step 3. (a) Compute the molar holdup A:‘). (b) Check for convergence. If converged, go to Step 5. (c) Set A(“)=A”~‘+(A!“)-_A(“-‘) / I J J 1.1 calStep 4. Perform the overall mass-balance culations to update A, and X, and go to Step 1. using
(45)
a,=O, a, = 1,
'dt
Step 1. Compute the bubble composition and temperature of the reboiler. x, = Y,.
form:
A.%=G(Y).
.
the stage flows Step 5. Calculate stage energy balances.
has the following
the
where
’
J, is the ith row of the Jacobian
matrix
of
G W. 500
=
450
E 4 12 0 400 .c &I .C :: ($
350 Relax * 1 3oc I
III1 01234
1 5
II
II 6
7
6
1IIlI 9
10
11
12
13
14
Iterations Fig. 2. Effect of relaxation factor on the calculation initial conditions.
of the
P. E. CUILLE and G. V. REKLAITIS
394 The above equations trivial equations,
can be combined
dyp
with the
dy,
dt=dt> for the differential
variables
to obtain
the system
G*=J*dy. dt
where J* is a matrix whose ith row consists of the coordinate vector eT if the ith equation is differential and of J,, if it is algebraic; and where, g: = dyi’)/dt if the ith equation is differential and g: = 0, if thejth equation is algebraic. This relation is a linear algebraic system, in the vector of unknowns dY/dt. If the J* matrix is nonsingular, the solution of this system will give all the initial derivatives. In the case of the system of equations (10)(29), the matrix J* is singular (the system has an infinity of solutions, including dY/dt). One method for solving this problem is to estimate a subset of derivatives in order to reduce the dimension of the system and to obtain a resulting nonsingular system. For instance, it is sometimes possible to assume that some derivatives are equal to zero at the initial time. However, this method is ad hoc and can not be guaranteed to work in all cases. Moreover, this kind of approximation is not possible with chemical reaction since with chemical reactions the initial derivatives will, in general, not be zero. Therefore, the following approach is preferable. The mathematical cause of this numerical problem is the discontinuity of the time derivatives of the solution functions at the starting point. This discontinuity has no physical reality since the start of production corresponds to the opening of a valve, which physically requires a small but nonzero duration. We can simulate the valve opening, by using a continuous function for the corresponding flow rate which increases the flow from zero to a finite positive value, in a small but nonzero time interval. The time derivatives of this function will be equal to zero at the beginning and at the end of this interval. For instance, if the column is to be simulated with constant production rate D = Do, then over the arbitrary time interval from, say, t = - 1 s to t = 0 the variable D in equation (18) can be replaced with the cubic polynomial function: D = D,( -2t3
- 3t2 + 1);
at
t=
-1,
D=O
and
g=O dt
and
t =O,
D =D,
and
dD -=O. dt
At the beginning of this time interval, i.e. at t = - 1, the column is operating at total reflux: all the time derivatives are equal to zero and are continuous functions.
Over the time interval - 1 to 0, the model equations are integrated with the cubic forcing function for D inserted into the model. At time t = 0 appropriate values of all derivatives of algebraic variables will have been calculated through the normal execution of the integrator and, thus, the integration can continue without interruption. In this way both the integration at t = - 1 s (total reflux condition) and at t = 0 (beginning of batch production) can be initiated without any numerical difficulties and without explicit estimates of the derivatives of algebraic variables. It should further be noted that this device can also be employed if other column specifications are of interest. For instance, if the column heating rate is specified as fixed and the distillate rate is left variable, then the reflux ratio variable can be subjected to the cubic forcing function. Specifically, at t = - 1, l/R = 0 since the column is at total reflux and at t = 0, l/R will take on the specified value. The transition function for l/R will thus correspond to the opening of a valve activated by a proportional controller. The only minor error introduced using this transition device is the fact that a distillate quantity equal to the integral of the distillate production curve over the time interval - 1 to 0 will have been produced prior to time 0. However, since the transition time interval can always be chosen to be small compared to the total batch distillation time, this error can be made as small as desired. (In our examples, we used 1 s for a 9000 s operating time.) EXAMPLES
OF SIMULATION
A general program simulating a batch rectification column has been developed, using the model and the numerical methods presented above. The program runs on a Control Data “Cyber 205” supercomputer and is written in FORTRAN 77 extended for vector calculations. This program can be easily adapted to any physical properties package, using small interface subroutines. The regular version uses the physical properties package of the simulator FLOWTRAN [8] through an interface which converts American units to SI units and adds heats of formation to the predicted enthalpies when reacting systems are indicated. Several examples of rectification with and without chemical reactions have been simulated, and will be reported next. In the first two examples, the options of the FLOWTRAN physical properties package used are the following: Antoine vapor pressure; Redlich-Kwong vapor fugacity; corrected liquid fugacity; ideal solution activity coefficients. Example
I
We have simulated successfully the rectifications performed by Domenech and Enjalbert [9] involving
Simulation of a staged batch distillation unit
Production
(mol)
Fig. 3. Simulation Example I. Experimental points and simulated curve for R = 3 and x,,, = 0.62 for C,H,,.
9000
Time
(s)
Fig. 6. Simulation Example time and plates for
1. Mole fraction
of C,H,,
vs
R = 2 and x,,, = 0.55.
the separation of cyclohexane and toluene. The experimental equipment used was a perforated-plate column, with 4-10 plates with a 60 1. reboiler heated with a heat transfer rate of 3 kW; this column operates under atmospheric pressure with a reflux ratio between 3 and 10 (see Ref. [9] for further details of the operating conditions). Figures 3-5 show the results of the simulations (curves) and the experimental data with reflux ratios of 3, 4 and 6 for different initial composition (respectively: molar fraction in cyclohexane of 0.62,0.55 and 0.42). Figures 68 show the variations of the
1
9000
o-
09- tb
Time
uou
080.7
6 06-
(s)
Fig. 7. Simulation Example time and plates for
1. Mole fraction
of C,H,,
vs
R = 4 and x,,, = 0.55.
-
E 058 G
04-
x” 0302Ol00
I 0
100
50
Production
150
200
(mol)
Fig. 4. Simulation Example 1. Experimental points and simulated curve for R = 4 and x,,, = 0.55 for C,H,,.
Time
(sl
Fig. 8. Simulation Example time and plates for
ooi-----k
Production
I 100
I 150
(mol)
Fig. 5. Simulation Example 1. Experimental points and simulated curve for R ==6 and x,,~ = 0.42 for C,H,,.
1. Mole fraction
of C,H,,
vs
R = 6 and x,,, = 0.55.
composition profile of cyclohexane in the column during the simulation, on a given example with reflux ratios of 2,4 and 6 (with the same initial cyclohexane molar fraction: 0.55). The flat portion of the surfaces drawn in Figs 6-8 corresponds to a zone in which the column performs a good separation. It is evident that the size of this flat zone increases with the reflux ratio, but this increase is achieved at the cost of higher utilities requirements.
P. E. CUILLE and G. V. REKLAITIS
396
Example 2 We also tested simulation cases with chemical reactions: for example the esterification of acetic acid with I-propanol which also involves the undesirable hydrolysis reaction: CH,-COOH
+ CH,--CH+ZH,OH
+ CH,-CQO--CH2-CH,-CH,
+ H,O
(46)
and CH,-COO-CH,-CH,--CH, + CH,-COOH
+ H,O
+ CH,-CH,-CH,OH.
(47)
We calculate the initial state at total reflux with a mixture of acetic acid and water (therefore without chemical reaction). Then we start the simulation in a fed batch mode: we feed the reboiler with a mixture of propanol and water and simulate the rectification with the production of propyl acetate. Figure 9 shows the molar fraction of ester vs the plate numbers and time. At t = 0, there is no ester in the column. In the reboiler [Plate 6) the reactions start immediately, but at the top of the column (condenser = Plate 0) the start of reaction is delayed. In the top plates, where the molar fraction of water is important, we can see that the fraction of ester vs time reaches a maximum, and then decrease, because of the hydrolysis reaction. Figure 10 shows the variation of the molar fraction of water. At t = 0, the column contains only acid and water which is the most volatile compound and the fraction of water is maximal at the top of the column. Just after the start of the rectification, because of the propagation of alcohol in the column, the fraction of water on the top plates falls down very quickly. Later, the fraction of water becomes maximum on Plate 2 because the alcohol is now the most volatile compound.
gooo sw Fig. 10. Simulation
Example 2. Mole fraction time and plates.
The initial parameters
of water
for this simulation
vs
are:
stages: reboiler, 5 plates, condenser; holdup of the condenser and of each plates = 0.25 dm3; Murphree efficiency of the plates = 70%; reflux ratio = 3; top pressure = 101,300 Pa, no pressure drop between plates; initial charge = 200 mol (acid 90%, water 10% mole fractions); composition of the feed = alcohol 90%, water 10% (mole fractions); feed rate = 0.02 mol SK’; esterification rate constant = 2 x hydrolysis rate constant (independent of the temperature). Example 3
p 5
Fig. 9. Simulation
Example 2. Mole fraction and plates.
of ester vs time
As a further illustration, we present a case involving the linkage of the simulation model to the more extensive thermodynamic properties package described by LeLann et al. [lo, 111. The simulation was executed on a Bull SPS9 computer operating under UNIX. The example involves the separation of the ethanol-water azeotrope from an initial mixture containing 10% (mol) ethanol in water. For this application, the liquid-phase behavior is modeled using the UNIQUAC method and the vapor-phase behavior using the Peng-Robinson equation of state. The results (not shown here) confirm that for operation of the column at 1 atm, one can, of course, obtain reasonably accurate results using the ideal gas equation of state. However, near the azeotropic condition, Peng-Robinson is clearly more accurate. Figure 11 shows the liquid composition profile for ethanol at the initiation of rectification and after 30 min of production in a 20-plate column. As might be expected, the initial flat profile tending towards the
391
Simulation of a staged batch distillation unit 1
.o
r
09 0.8 07 06 '05 0.4 0.3 02 0.1 0.0 0
1 2 3 4 5 6
7 8 9101112131415161718192021 Plates
Fig. 11. Simulation Example 3. Mole fraction of ethanol vs plates.
azeotropic composition becomes steeper over time as the alcohol-rich distillate is removed. USE OF THE “CYBER
205” SUPERCOMPUTER
The main advantage of this type of supercomputer, compared with the classical “Von Neuman” machine, is the capability to perform parallel calculations when executing operations on vectors. On the Cyber 205, two ways of using this capability are available:
I. The automatic
“vectorization” of the do-loops: if a do-loop satisfies certain vectorization constraints, the Fortran compiler can translate it into a machine vector instruction. If those vectorization constraints are satisfied, programs written in standard Fortran 77 (or Fortran 66) can use the capabilities of the vector processor. 2. The explicit vector programming: using some extensions of the Fortran language (specific to the Cyber 205) it is possible to program directly in a vectorized form. The integration package used in the simulation was written in standard Fortran 77, however, it is not well-adapted to the vector capabilities of the machine. Since only a few of the do-loops are vectorizable, the use of the automatic vectorization of the loops on a given example of simulation (the cyclohexane-toluene separation) saves only 3% of the CPU time compared with a simulation without this vectorization. In order to enhance the use of the vector capabilities, we rewrote the solver used by the integrator at the corrector step. The original one used a Gaussain elimination method (GEM) and made very poor use of the vector processor. The new one uses a Jordan method which is numerically less efficient on a classical computer, but is better adapted to be programmed in a vectorial form. For instance, on a separate example involving a system of 50 algebraic linear equations, this routine is more than four times faster than that based on a GEM. Using this routine in the simulation, we can save more than 40% of the CPU time on the cyclohexane-toluene example.
The actual version of the simulation program contains 2081 lines of Fortran code (without the integrator and the FLOWTRAN physical properties package). The CPU time needed for a simulation depends on the characteristics of the simulated problem. For the same column (5 plates) the program need only 3 s of CPU for simulating 9000 s of operation for the conditions of Example 1 (two compounds, no chemical reaction), and 73 s of CPU for the conditions of Example 2 (four compounds, two chemical reactions) also for 9000 s of operation. A better use of the vector processor is limited by a major bottleneck: the calculation of the physical properties of the gas and liquid mixtures used in the rectification. There is no physical properties package written in a vectorized form available today. If the use of computers with vector capabilities become more widespread, it will be necessary to develop new packages able to calculate not only one physical property at a time, but also a series of values of a given property for a set of different mixtures. For example: to calculate simultaneously all the liquid enthalpies for all the liquid mixtures present in the rectification column at a given time. Acknowledgemenrs-This work was carried out with the financial support of Rhone Poulenc, S.A., which is gratefully acknowledged. The work was completed at Purdue University while the first author was on leave as a visiting scholar and is part of the research presented in his Docteur Thesis.
NOMENCLATURE
A D e F H, H,
= = = = = =
K= k = L = M = n = n, = neq = P = $ 1
Material holdup (mol) Overhead product flow rate (mol s-‘) Murphree efficiency Feed rate Molar enthalpy of the liquid phase Molar enthalpy of the vapor phase Vector of k-values Element of the vector K Liquid flow rate (mol s-‘) Molar mass Number of plates Number of components Number of equations Pressure (Pa) ;;Fu;frit;tt transfer (W)
R = Vector of chemical reaction rates s = Iteration counter T = Thermodynamic temperature (K) t = Time(s) V = Vapor flow rate (mol s-l) X = Vector of mole fractions (liquid) x = Element of the vector X Y = Vector of mole fractions (vapor) y = Element of the vector Y An = ZR, r = Volume (m’) p = Density (kg mm3) 1 = Relaxation factor Subscripts 0 = Condenser 1 = Plate 1 (top) j = Plate j
P. E. CUILLE and G. V. REKLAITIS
398 n b i ini f
= = = = =
Plate n (bottom) Reboiler Component i Initial Feed
REFERENCES 1. E. L. Meadow, Chem. Engng Prog. Symp. Ser. 46(59), 48 (1963). 2. G. P. Distefano, AIChE JI 14, 190 (1968). 3. C. W. Gear. Communs Ass. comout. Mach. 14. 176 (1971). 4. J. F. Boston, H. I. Britt, S. Jirapongphan and V. B. Shah, An Advanced System for the Simulation of Batch
Distillation Operation, Vol. 11, pp. 203-237. Foundation of Computer Aided Chemical Process Design, Engineering Foundation, New York (1980). G. D. Byrne, A. C. Hindmarsh and K. C. Jackson, _J. Comput. ihem. Engng 1, 133 (1977). 6. S. E. Gallun and C. D. Holland, Comput. them. Engng 6, 231 (1982). Ass. comput. Mach. Signum News1 15, 7. A. C. Hindmarsh, 10 (1980). 8. J. D. Seader, W. D. Seider and A. C. Pauls, FLOWTRAN Simulation, an Introduction, 2nd edn. CACHE, Ann Arbor, Mich. (1977). 9. S. Domenech and M. Enjalbert, Chem. Engng Sci. 29, 1519 (1974). Docteur Ingenieur Thesis, INP 10. J. M. LeLann, Toulouse, France (1983). 11. J. M. LeLann, X. Joulia and B. Koehret, Entropic 123, 29 (1985).