Ind. Eng. Chem. Res. 2008, 47, 813-824
813
Startup Operation of Middle-Vessel Batch Distillation Column: Modeling and Simulation Sven Gruetzmann* and Georg Fieg Hamburg UniVersity of Technology, Institute of Process and Plant Engineering, Schwarzenbergstr 95c, 21073 Hamburg, Germany
The objective of this article is to extensively analyze and discuss the startup of a middle-vessel batch distillation (MVBD) column by means of dynamic modeling and simulation. A rigorous model describing the entire process of an MVBD column including startup from ambient conditions is presented. The model considers the heating of the distillation column and peripheral equipment, as well as the formation and propagation of hydrodynamic and thermodynamic profiles. The physical validity is shown by comparison of simulation studies with experimental results. Sensitivity analyses are used to study startup operation and promote a deeper understanding of the process. General guidelines are provided to standardize and improve the process. The theoretical examinations offer the basis for the development of improved process control strategies and recipebased automation concepts, as well as for dynamic optimization studies. Introduction Batch distillation is a standard unit operation in the production of specialty and fine chemicals, as well as pharmaceuticals. These industries are characterized by the fact that their product portfolio is composed of high-purity and high-value-added products. Therefore, it is reasonable that batch distillation has been under renewed investigation in the past two decades. Since the beginning of the 1990s, investigations have focused on novel batch distillation columns, namely, middle-vessel1-4 and multivessel5-8 columns. A representation of an MVBD column is shown in Figure 1. Both the location of the feed and the process control strategy are free to choose. Actually, the process was designated to separate light- and heavy-boiling components from a middle boiler by charging the feed to the middle vessel and removing the impurities at the top and bottom of the column.1 Furlonge et al.7 found that little economic incentive exists to optimize the initial feed location and recommended supplying the feed to the reboiler. This is reasonable from a practical point of view as well. A favorable process control strategy for high-purity products from multicomponent mixtures has been discussed by various authors.2,5,9 As products accumulate in the side vessels and no stream is removed from the column, the process approaches a steady state. In this case, no product streams, denoted LD and LB, exist (Figure 1). If the products are within specifications, they are discharged. If not, set points can be changed to achieve the desired specifications as long as the separation performance is sufficient. The column is operated with infinite reflux and infinite reboil ratio and is called closed operation or total reflux operation. Because this approach exhibits a number of significant advantages (operation with maximum separation performance, simple process control, no off-cuts), this article considers the more attractive closed operation mode. Most publications dealing with the modeling and simulation of distillation columns neglect the startup period and its influence on the process. Experimental studies by Gruetzmann et al.10 recently indicated that the deliberate switching of discrete decision variables (e.g., reflux flow) can lead to a shorter * To whom correspondence should be addressed. Tel.: +49 40 42878-3241. Fax: +49 40 42878-2992. E-mail:
[email protected].
process. In terms of optimization, a focus on the first period (see Figure 1), in which hydrodynamic and thermodynamic profiles are formulated, is reasonable. The startup of continuous distillation has been covered by several works. 11-14 An analysis of the startup period of conventional batch distillation columns with product removal has been considered in only a few cases. For example, Wang et al.15 and Elgue et al.16 presented detailed models simulating the startup of batch distillation. However, the total reflux operation of novel batch distillation columns has not been considered. Gruetzmann et al.17 briefly presented an appropriate model for investigating the startup of an MVBD column for the first time, but without going into detail. This article extends the previously mentioned work by using the model to analyze the startup period of an MVBD column. This article is divided into two major parts. In the first part, we present a mathematical model capable of handling physical phenomena occurring during startup. In addition, we address the solution procedure and solver handling. We show that the model can accurately describe the physical phenomena occurring on a theoretical stage. Then, we validate the model by comparing experiments performed on a laboratory-scale glass column with simulation results. The second part considers the influence of the startup period on the process. Simulation studies are used to show the difference between the developed model and the common “simplified” approach with initially filled trays at the boiling temperature. Sensitivity analyses were performed to illustrate the effect of chosen parameters, specifically, the overall mass of the column, heat losses, internal column holdup, amount of feed, and changes in feed composition. This article concludes by providing the first practical heuristics for improving and generalizing the process. This elaborated methodology provides a strong basis for a deeper understanding of the process. The information will be used to develop advanced process control strategies and recipe-based automation concepts and to carry out dynamic optimization studies including startup in future works. Modeling Approach The inherent change of state variables in batch distillation is typically highly nonlinear, particularly during startup. The corresponding difficulties in terms of initialization of variables
10.1021/ie070667v CCC: $40.75 © 2008 American Chemical Society Published on Web 01/09/2008
814
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008
Figure 1. Process scheme and diagram of a middle-vessel batch distillation column.
Figure 2. Representation of tray j.
and high stiffness of the resulting system of differential/algebraic equations (DAE system) lead to strong requirements of numerical solvers. This is the reason that mostly shortcut methods or simplified “rigorous” models have been used. Common simplifications are, for example, the assumption of constant relative volatilities, constant molar liquid holdups on all theoretical stages, or constant heats of vaporization. The initialization is usually performed assuming filled trays with the feed at the boiling temperature. The advantage of these simplified models is, of course, that much information on the column dynamics can be obtained with relatively little effort. However, in the case of batch distillation, the initial distribution of the feed during startup is important for the subsequent course of optimal trajectories of the manipulated variables. Because the formation of the composition profiles depends on both the overall heat transfer and the heating of the column, it is recommended that these effects be investigated. The model developed here is intended to be used for sensitivity analyses, the development of process control strategies, and dynamic optimization studies. Furthermore, we intend to use the model for the efficient planning of experimental studies on a pilot-plant distillation column. The latter will be operated under vacuum and therefore contains several sections of structured packings that can often be found in the fine- and specialty-chemical industries as well. The height equivalent to a theoretical plate (HETP) value has been used to set up a tray-
to-tray equilibrium model in the modeling and simulation software Aspen Custom Modeler. In this section, only the mathematical model for a theoretical stage is presented. Some modifications are made for the reboiler, accumulators and distributors, condenser, and product tanks. To describe the startup behavior accurately, the following features were included in the model: variable relative volatilities, ideal and nonideal vapor-liquid equilibria, variable molar vapor flow, variable molar liquid holdup on the trays, nonadiabatic operation, and variable pressure drop. Thus, the well-known MESH equations were extended by the heat accumulation of the process equipment and complemented by algebraic relationships calculating the transient tray hydraulics, the maximum molar liquid holdup, the pressure drop, and possible heat losses. The necessary physical and thermodynamic properties of the components and vapor-liquid equilibrium data were provided by internal procedures from Aspen Properties. Figure 2 shows a representative tray with the main parameters. Material and Energy Balances. Index i denotes the components, and index j denotes the theoretical stages. The numbering is from the top to the bottom. The vapor holdup on a tray can be neglected because the vapor and liquid molar densities differ by several orders of magnitudes in vacuum applications. Thus, the overall material balance and component balances can be written as
dHUj ) Lj-1 + Vj+1 - Lj - Vj dt
(1)
d(HUjxi,j) ) Lj-1xi,j-1 + Vj+1yi,j+1 - Ljxi,j - Vjyi,j dt
(2)
The energy balance includes an accumulation term for the equipment and takes heat losses into consideration.
d(HUjhj,liq + Hj,internals + Hj,shell) ) Lj-1hj-1,liq + dt Vj+1hj-1,vap - Ljhi,liq - Vjhj,vap - Qj,loss (3) Hj,internals + Hj,shell ) (mcpT)j,internals + (mcpT)j,shell
(4)
The heat capacities of the column internals and the shell are considered through the quantity mcp with indices “internals” and “shell”. Heat losses are determined from eq 5. The overall heat-transfer coefficient k is calculated in advance using common correlations. The surface area is given by the column design.
Qj,loss ) kA(Tj - Tambient)
(5)
2
[ (
a2/3 1 - Fvap (ν u )1/3 1 C 2 3 dp liq liq ∆p ) 2 L 1 - Fvap a1/3 2/3 -5 ψ 3 1 - C1 u dp liq ψ
)
∑i xi,j ) 1
(6)
∑i yi,j ) 1
(7)
Vapor-Liquid Equilibrium. If both vapor and liquid streams exist, vapor-liquid equilibrium is calculated from
yi,j ) Kixi,j
(8)
where K denotes the equilibrium constant and is called from Aspen Properties depending on the chosen physical equilibrium method. Because the column hydrodynamics vary frequently during startup, the actual HETP value differs from its steadystate value. Although changing the structure of the model and re-initializing during the simulation is not possible, Kruse18 introduced a Murphree efficiency coefficient that can take various values to overcome this problem. However, Flender19 showed that the impact on the results is very small. Within our modeling approach, the transient deviation from equilibrium is neglected. Tray Hydraulics. The first vapor arises if the temperature in the reboiler approaches boiling temperature or, in other words, if the vapor pressure of the mixture equals the pressure in the reboiler. This is also true for each equilibrium stage. However, in the model, the vapor pressure is called from an Aspen Properties subroutine and cannot be set equal to the tray pressure because it causes an index to be greater than 1. Therefore, we defined a dummy variable for the vapor pressure and calculated the vapor stream from the equation
pVL,j - pj pj
]
-5
∀ Reliq < 2 ∀ Reliq g 2 (11)
with
ψ ) K1RevapK2
Summation Equations.
Vj ) C
{
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 815
Revap )
uvapFvapdp (1 - )ηvap
Reliq )
uliqFliqdp ηliq
(12) (13)
(14)
The influence of the calculation method was investigated by comparison of two simulation runs. In the first simulation, the pressure drop was calculated in every time step using eqs 1114. In the second one, the total pressure drop was taken from the previous run, and a linear pressure profile was assumed. Figure 3a shows the propagation of the pressure profile. The total pressure drops in the two cases were equal, as indicated by the dotted line. Figure 3b shows the temperature course on selected trays. The total number of trays was 16.
(9)
If the vapor pressure of the mixture, pVL,j, is larger than the pressure on the tray, pj, a vapor outlet stream, Vj, is computed. The constant is adjusted to reduce the difference between that pVL,j and pj. The tray pressure is determined using the equation
pj ) pj-1 + ∆pj ) f(geometry, physical properties, tray hydraulics) (10) The pressure drop, ∆pj, is an important parameter in vacuum distillation, because the overall pressure drop is of the same order of magnitude as the tray pressure. It is defined as the difference between the pressures on two adjacent trays. Special correlations from the open literature and from vendors have been used to calculate the pressure drop for different types of structured packings.20,21 As an example, the following equations were used to calculate the pressure drop by Mackowiak.21 We refer to the Nomenclature section for an explanation of the symbols.
Figure 3. Influence of the pressure calculation method on the process.
816
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008
Figure 4. Comparison of maximum holdup calculations on a representative tray.
It can be concluded that the influence of the pressure drop calculation method is rather small. Therefore, we recommend using correlations to determine an accurate steady-state pressure drop in advance and employing a constant value in simulations to increase the robustness of the model. The liquid outlet stream is calculated in a similar way (eq 15). The parameter τliq represents the hydraulic time constant and has to be adjusted depending on the column that is investigated.
Lj )
HUj - HUj,max τliq
(15)
The maximum liquid holdup on a tray can also be calculated by different correlations.20,21 It depends on the geometry, the physical properties of the liquid phase, and the vapor and liquid loads. It is the sum of the static and dynamic holdups.
HUj,max ) HUj,stat + HUj,dyn ) f(geometry, physical properties, tray hydraulics) (16) We carried out simulation studies to investigate the influence of the dynamic calculation of the maximum holdup in accordance with Engel et al.20
[
HU ) HUstat + HUdyn0 1 + 36
( )] ∆p/L Fliqg
2
HUstat ) 0.033 exp(-0.22Bo) Bo )
HUdyn0 ) 3.6
g1/2
(18)
gFliq
(19)
σliqa2
( )( )( ) uliqa
0.66
ηliqa3/2 Fliqg1/2
(17)
0.25
σliqa2 Fliqg
0.1
(20)
The results are illustrated in Figure 4. The solid lines represent the case in which the maximum holdup was employed in the
dynamic simulation using eqs 17-20. The dashed lines correspond to a simulation run in which the maximum steady-state holdup was calculated in advance using the same correlation. The first vapor enters the tray after approximately 150 s. Because of condensation, the liquid holdup increases. If the holdup is calculated dynamically, the maximum liquid holdup equals the static holdup until a liquid stream enters the tray after approximately 500 s. The calculated liquid and vapor outlet streams are slightly delayed, but no influence on the process can be observed. It can be clearly seen that there is no need to include the holdup calculation within the dynamic simulation. It is recommended that one proceed in the same way as in the pressure drop calculation, that is, determine the maximum liquid holdup in advance and set a constant value during simulations. The overall process model consists of similar sets of equations to describe the reboiler, the condenser, the distillate receiver, and the middle vessel. Additional control models are provided by the software package and can be easily implemented to the model. Solution Procedure An appropriate selection of the numerical solver depends on the characteristics of the DAE system. The main criteria are the stiffness and the index of the system, as well as the occurrence of discontinuities. A system is called stiff if it contains time constants that vary by several orders of magnitude. This is true for the system investigated in this work. As a result, only implicit solvers can be considered. The index of a DAE system is the number of times one needs to differentiate the algebraic equations to convert the system to an ordinary differential equation (ODE) system.22,23 Solvers provided by Aspen Custom Modeler are limited to simulations of DAE systems with index 1. Mathematical modeling of physical systems, however, can result in index larger than 1, for which no general-purpose solvers are available.24,25 Two general techniques can be used to reduce the index of a system. The first is the transformation of a DAE system into an ODE system by differentiating the appropriate equations. An analysis of the
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 817
Figure 5. Behavior of tray j during startup.
present DAE system shows that the equations that cause the index problem are internal subroutines. Thus, this solution method is not applicable in our case. The second approach is the introduction of dummy variables and additional equations. The successful application of such a procedure is described in the Modeling Approach section. The third challenge that the solver must handle is discontinuities that always occur when an outlet stream is computed for the first time, that is, when the value changes from zero to nonzero. Such discontinuities are typically indicated by IF/THEN or STATE/SWITCH statements within the programming code (see also Wang et al.15). The following three implicit solvers are provided by Aspen Custom Modeler to integrate the DAE system: Gear algorithm, implicit Euler, and VSIE 11.1 (variable-step implicit Euler). All of these methods were successfully applied. In the following studies, the implicit Euler method with a variable step size was used for the integration step. A minimum step size of 1 × 10-5 h was chosen to guarantee sufficient accuracy. The Newton method was employed to solve the nonlinear algebraic equation system. The Jacobian matrix was calculated at every iteration to ensure convergence because of the stiffness and the occurrence of discontinuities. Simple linear equations were solved by the solver MA48 using default settings.26
column is 20 °C. No stream enters or leaves the empty tray j. During period 1, a vapor stream Vj+1 enters the empty tray. With the vapor condensing at the column wall and internals, the tray is filled to its maximum liquid holdup, HUmax. If the liquid holdup HUj exceeds the maximum holdup, a liquid stream Lj flows down the column during period 2. The rising vapor stream still heats the tray until the liquid reaches boiling temperature. At this stage, a vapor-liquid equilibrium is calculated. If the vapor pressure pVL,j is larger than the pressure on the tray pj, a vapor outlet stream Vj is calculated (period 3). Period 4 indicates that liquid from the upper tray flows downward (Lj-1). Note that periods 2-4 can appear in a different order depending on the maximum liquid holdup, the strength of the condensing effect, and the boilup rate. Assuming a high reboiler duty and a nearly adiabatic process, the following observation can be made: The temperature on a tray increases until a vapor stream is computed (period 3). Because only a small amount of liquid accumulates on the tray, the tray is filled by the internal reflux from the top of the column (period 4). Afterward, a liquid outlet stream is generated (period 2). We can summarize that the simulation results clearly show the physical consistency of the model. Model Validation
Physical Validity The objective of this section is to verify that the dynamic model represents the entire process, beginning from a column at ambient conditions. The term “ambient conditions”, therefore, means that the equipment is cold, the feed is entirely distributed in the product tanks, and the trays are empty. The column is heated by the rising vapor. The trays are filled by condensation and reflux flow. A consistency check shows that the physical phenomena are described correctly, despite the fact that the presented model is a discrete representation (theoretical stage) of a packed column. Hence, the typical courses of the temperature, the internal molar flows, and the liquid holdup of a representative tray j during startup are shown in Figure 5. At the initial state (period 0), the column is at ambient conditions, that is, the temperature of both the feed and the
In this section, we describe the application of the dynamic model to a laboratory-scale glass column the objective of which was to separate a ternary mixture of fatty alcohols into three pure products. The experiments were carried out in cooperation with an industrial partner (Cognis Deutschland GmbH & Co. KG, Du¨sseldorf, Germany). The main parameters are summarized in Table 1. Both the liquid holdup on a tray and the pressure drop were calculated in advance using the procedure described above. Information provided by the vendor was used to estimate the HETP value. The overall heat-transfer coefficient (kA ) 0.02 W m-1 K-1) was determined using fundamental heat-transfer equations, and the column weight was estimated on the basis of the geometry. A comparison of the simulation results with the experimental results verifies the validity of the model. The MVBD column was operated with infinite reflux
818
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008
Figure 6. Comparison of simulation and experiment. Table 1. Main Parameters of the Laboratory-Scale Column and Simulation Setup parameter
simulation
experiment
column diameter type of packing N (upper section) N (lower section) feed mixture feed charge composition reboiler duty top pressure pressure drop column holdup
70 mm 6 5 hexanol/octanol/decanol 1665 kg 24.3/25/50.7 (wt %) 1000 W 0.2 kPa calculated calculated in advance
70 mm 2 × 1 m Sulzer BX 18.6 cm 15.5 cm hexanol/octanol/decanol 1665 kg 24.3/25/50.7 (wt %) 1000 W 0.2 kPa 2.5 mbar (overall) -
and reboil ratio. No products were removed until the end of the process. The feed was charged to the reboiler, and distillation was performed under vacuum. Figure 6 shows the temperature at different locations as a function of the relative process time, which is given by the actual time divided by the final process time. Tray 1 and the reboiler represent the top and bottom of the column, respectively. The temperature measurements used for process control were made on trays 4 and 10. The simulation results agree very well with the experiments despite the high level of discontinuities and nonlinearity. The deviations can be explained by some column data that could not be measured. For example, the overall column weight was estimated. Another important uncertainty is the choice of the inverse HETP value, so that the number of theoretical stages was set to 5 for the upper column section and 6 for the lower column section. This is the reason for the slight steady-state deviation of the temperature on tray 4. We can summarize that the results demonstrate the validity of the dynamic model. Simulation Studies Every batch distillation is characterized by frequent startup and shutdown operations, although, in most cases found in the open literature, the startup of simple or complex batch distillation column has not been considered.5-9 This is because startup usually comprises only a small part of the entire process. In the second part of this article, we show why it is beneficial to examine this short period. Our objective is to analyze the process of an MVBD column and try to make the first general suggestions for process control based on our analysis. In future works, improved process control strategies and recipe-based automation concepts will be developed on the basis of the detailed theoretical examinations. Because transfer of process control strategies developed for small columns might result in poor performances in industrial applications, we consider most
Figure 7. Temperature control in closed operation mode. Table 2. Column Dimensions and Other Parameters parameter
value
column diameter N (upper section) N (lower section) feed mixture composition
68 mm 8 8 hexanol/octanol/decanol 33.3/33.3/33.4 (wt %)
types of columns by varying relevant parameters. The main parameters influencing startup are heat transfer to the surroundings, heat transfer to the equipment, internal holdup, feed charge, and feed composition. Temperature control using PI controllers was applied in all studies, that is, the process was characterized by an upper control loop and a lower control loop (Figure 7). The two temperature set points were determined by the arithmetic mean value of the boiling temperatures of the pure components separated in the upper and lower column sections. The column dimensions and other parameters used for simulation can be found in Table 2. Product specifications were subjected to tight control, that is, wi > 0.99 kg kg-1. Note that, for the second part of this article, the column dimensions were taken from the pilot-plant column at the Institute of Process and Plant Engineering (Hamburg University of Technology) and, therefore, differ from the data used in the Model Validation section. Several cases were investigated, as summarized in Table 3. In all cases analyzed, the reboiler duty was adjusted after startup to guarantee similar maximum vapor loads at the top of the column. Cases A-C refer to the startup model presented in this work. Both the column and the feed mixture are at ambient conditions. In cases D and E, initialization was adjusted so that startup was neglected. Then, the following conditions were assumed: trays filled with liquid at boiling temperature, feed at boiling temperature, no pressure drop, and adiabatic operation. In case D, the initial composition on the trays was that of the feed, whereas in case E, the liquid on the trays contained pure low boiler. In the following discussion, these cases are called “simplified initialization”. Cases A-C can be distinguished by
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 819
Figure 8. Temperature profiles during the operation of an MVBD column. Table 3. Main Parameters for Cases Considered A
B
parameter
units
1
2
3
4
HUF TF Qreb (initial) Qreb (steady-state) p ∆p (steady-state) q1 q2 HUmax HUj (initial) Tj (initial) xj (initial)
kg °C W W kPa kPa W/W K/s mL mL °C
7.5 20 2000 2000 1.5 1.2 0.6 0.03 177 1 × 10-6 20 xFeed
25 20 2000 2000 1.5 1.2 0.6 0.03 177 1 × 10-6 20 xFeed
75 20 2000 2000 1.5 1.2 0.6 0.03 177 1 × 10-6 20 xFeed
7.5 20 2000 2000 1.5 1.2 0.6 0.03 10 1 × 10-6 20 xFeed
C
7.5 20 2000 800 1.5 0.6 0 0.03 177 1 × 10-6 20 xFeed
1
2
3
4.8 20 2000 800 1.5 0.6 0 0.1 177 1 × 10-6 20 xFeed
25 20 2000 800 1.5 0.6 0 0.1 177 1 × 10-6 20 xFeed
75 20 2000 800 1.5 0.6 0 0.1 177 1 × 10-6 20 xFeed
D
E
4.8 75 800 800 1.5 0 0 177 177 75 xFeed
7.5 61.5 800 800 1.5 0 0 177 177 61.5 xC6 ) 1
Table 4. Simulation Results A parameter startup temperature control (distillate receiver) temperature control (middle vessel) termination criterion 1: all components within specifications termination criterion 2: (99% of steady-state value) end of product accumulation process duration limiting component (termination criterion 1)
1
2
3
4
t1 t2 t3 tt1
h h h h
0.68 1.32 1.43 1.43
1.10 2.61 2.94 4.47
2.31 5.85 7.49 >15
0.65 1.32 1.46 2.04
tt2
h
1.61
4.34
11.75
t4
h
1.61 C10
4.47 C10
>15 C8
the overall heat-transfer ratio q1 (eq 21) and the heating rate q2 (eq 22).
q1 ) q2 )
∑Qloss Qreb
(W W-1)
Qreb mequipmentcp,equipment
B
units
(K s-1)
(21)
(22)
Note that q1 and q2 are, in fact, related to the initial reboiler duty. The individual cases are characterized as follows:
C
D
E
1.71 3.20 10.32 32.01
0.00 0.60 0.81 2.34
0.00 0.89 0.92 1.04
4.64
28.38
2.30
1.04
9.41 C8
32.01 C8
2.34 C8
1.04 C10
1
2
3
0.57 1.15 1.30 1.39
0.28 0.59 0.96 1.83
0.66 1.65 3.39 9.41
1.86
1.77
1.46
2.04 C8
1.77 C8
1.83 C8
Case A. Both heat losses and heat accumulation were taken into consideration. This case might occur for laboratory-scale columns; for example, Kruse18 mentioned heat losses of more than 50% of the reboiler duty for a laboratory-scale glass column. The value of the heating rate is rather small and typical for small columns made of steel or double-walled glass columns. Case B. Heat accumulation to the surroundings was considered, but not heat transfer. This case is typical for well-insulated pilot-plant columns, in which a nearly adiabatic process can be assumed.
820
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008
Case C. Heat losses were neglected, and the heat accumulation effect was low. This case is representative for well-insulated industrial-scale columns or well-insulated “lightweight” laboratory-scale columns. The heating rate q2 is increased for industrial-scale columns because scale-up of the reboiler duty and the overall column weight is usually nonlinear. Results in terms of significant times are summarized in Table 4. The time needed to heat the column is denoted as t1. The control loops in the upper and lower column sections are activated at t2 and t3, respectively. In practice, the process can be terminated if both product specifications are met and holdup in the product vessels has achieved steady-state values. These time points are indicated by tt1 and tt2, respectively. The entire process duration is defined as t4. Additionally, Table 4 reports the limiting component of the process. Startup Operation The first simulation run (case A1) is used to briefly explain the general process behavior. Figure 8 shows the temperature profile at various times in the period between startup and steady state for a simulation in which heat losses and heat transfer to the column shell were considered. The process can be subdivided into three parts: startup, filling of the reflux drum, and profile propagation. Startup. At the beginning of the first part of the process, the column is at ambient conditions, that is, it is cold, and the trays are empty (Figure 8a). Rising vapor heats the column and reaches the condenser after t1 ) 0.68 h (41 min). No hydrodynamics is established so far. Filling of the Reflux Drum. As condensed liquid accumulates in the distillate receiver during the second part of the process, the middle boiler evaporates from the reboiler, resulting in a rise of the temperature profile (Figure 8b). Profile Propagation. The upper set point is reached after t2 ) 1.32 h (79 min). The reflux valve is opened, and liquid flows down the column to fill the middle vessel until the second control loop is activated (t3). Figure 8c clearly shows the separation effect and the detailed formation of developed temperature profiles. After 180 min, the temperature profile approaches a steady state. Termination criteria are fulfilled earlier, that is, tt1 ) 85 min, and tt2 ) 97 min. Although the first two parts of the process are usually short compared to the third, it is obvious that the first control actions at the end of the second part depend on the preceding progress. To further accentuate the need for an accurate model, we compare the startup model presented in this work with a common approach in which trays are initialized with maximum liquid holdup, feed composition, and liquid at the boiling temperature. For example, Furlonge et al. used a similar approach for dynamic optimization studies of a multivessel batch distillation.7 Parts a and b of Figure 9 represent to cases C1 and D, respectively, and show the key composition responses in the product vessels for the two approaches. Note that the origin indicates the time at which vapor first enters the condenser. Two effects can be pointed out when using a simplified initialization. First, the process termination criteria are fulfilled later even if the startup period is subtracted as in case C1 (2.31 h compared to 1.83 h). The reason for this difference is that impurities caused by the simplified initialization are trapped in the product vessels. Second, the trajectories in the two cases take completely different paths to steady state. One might argue that trays can be initialized with liquid rich in low boiler to start the process with a more realistic assumption
Figure 9. Key composition responses in the product vessels.
as was done by Skogestad et al.6 Therefore, we examined the case in which feed is charged to the reboiler and the trays contain only hexanol (case E). The results are illustrated in Figure 9c. Indeed, the composition responses in Figure 9a and c are quite similar in shape for hexanol and decanol. However, the amount of impurities in the middle vessel is rather small compared to startup model, which leads to a different conclusion with respect to the limiting key component. Moreover, the process is faster compared to case C1 (see Table 4). It can be summarized that the startup model describes the process more realistically. This is very important for drawing the right conclusions from simulation results. Thus, we believe that using the startup model for the development of improved process control strategies, recipe-based automation concepts, and dynamic optimization studies is required.
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 821
Figure 10. Impurities in the middle vessel for cases A1 and A4.
Influence of Internal Column Holdup and Amount of Feed The amount of feed and the internal column holdup can change dramatically if different types of columns are examined. The feed charge in laboratory-scale columns is usually small, whereas the internal holdup is relatively high. In contrast, industrial-scale columns often supply a feed charge that is much higher than the column holdup. We carried out two sensitivity analyses to examine the influence of these parameters. First, we kept the amount of feed constant and varied the internal holdup. The internal holdup is the sum of the liquid adhering to the structured packings and the liquid accumulated in the distributors. A change in this parameter can be justified by the use of different types of structured packings or internals. Two cases were compared. In the first case, the relative amount of liquid within the column was approximately 36 vol % of the feed charge (case A1). Then, the holdup was reduced to 6 vol % while the other parameters were held constant (case A4). The results are summarized in Table 4. The duration of startup and the times at which the closings of the two control loops were performed are nearly identical (t1-t3). One can therefore assume that the entire process runs similarly. However, the processes terminated after 1.61 and 2.04 h, respectively. Moreover, the limiting component switched from decanol to octanol. Figure 10 shows the mass fraction of impurities in the middle vessel during the process. One can see that the amounts of low boiler entering the middle vessel are the same in the two cases, and therefore, this parameter is not the reason for the difference in process duration. In fact, the limiting impurity in the middle vessel is the high boiler, as the composition trajectories show. Here, the larger internal holdup in case A1 serves as a buffer against impurities in the middle vessel. This example shows the importance of a startup model with regard to a deeper understanding of the process, especially if the objective is to develop generalized process control strategies. In the second sensitivity study, we increased the feed charge (7.5, 25, and 75 kg) while maintaining the absolute internal holdup (cases A1-A3 and C1-C3). As a selection only, the mass fraction of the key components in the product vessels are displayed as a function of time for cases A1-A3 (Figure 11ac). It is obvious that the process time increases with an increase in feed charge. More remarkable is the observation that the process is significantly faster in cases A2 and A3 than in the related cases C2 and C3 (see Table 4). As explained, the reboiler duty in cases C1-C3 was decreased after startup to give comparable vapor loads at the top of the column in all cases. However, the heat losses considered in cases A1-A3 lead to
Figure 11. Responses of the mass fractions of the key components in the product vessels.
larger internal streams that are responsible for the faster process. This observation again emphasizes the need for an accurate model as presented in this work. Influence of Heat Accumulation and Heat Losses In this section, additional simulation results are presented to analyze the process and show the influence of the overall heattransfer ratio, q1, and the heating rate, q2. Figure 12 displays the temperature profiles for cases A1, B, and C1 at t1-t4. Temperature sensors for control purposes should ideally be located in the middle of a column section in accordance with Wittgens et al.9 Because each column section was discretized into eight theoretical stages, we calculated the average temperatures between trays 4 and 5 and between trays 12 and 13. The
822
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008
Figure 12. Temperature profiles of cases A1, B, and C1.
analysis is divided into three parts corresponding to the three significant parts of the process. Startup. The first conclusion is that no fully developed profiles exist before liquid is sent back from the reflux drum. Most of the column is filled with light-boiling component. This fact confirms what we have already observed. The second conclusion is that the profiles differ depending on the varied parameters. If more heat is transferred to the surroundings, the internal liquid flow is higher, which causes a separation between the low boiler and the middle boiler. The consequence is that the middle boiler remains in the lowest part of the column indicated by the low temperature level. Filling of the Reflux Drum. The shapes of the profiles for cases A1 and B are similar in the second part but clearly differ from that for case C1. As explained earlier, the internal reflux due to the strong condensation effect leads to a separation of the light and middle boilers. The typical s-shape is already formed in the upper column section. Because the reboiler duty is adjusted afterward, it takes approximately the same time to fill the reflux drum in cases A1 (t3 - t2 ) 0.63 h) and B (0.58 h). However, in case C1, it takes only 0.3 h until the temperature in the middle of the upper section approaches the set point, which causes control action. The reason for this difference is the large amount of middle boiler evenly distributed and indicated by the low gradient of the profile in the upper column section. The fact that cases A1 and B show similar behaviors leads to the conclusion that the influence of the heat transfer to the equipment indicated by q2 is strong in this case. Profile Propagation. If reflux from the distillate receiver is started (t2), the connection to the middle vessel is opened, and liquid from the upper column section accumulates in the middle vessel. The later temperature control is activated, the more impurities can accumulate in the middle vessel. The amount of low boiler trapped in the middle vessel depends on the parameters discussed above; for example, the amount is higher
in case C1 than in case A1. Therefore, the process needs more time to achieve steady state; that is, the third part (t4-t2) of the process is longer (1.24 h compared to 0.29 h). From a comparison of case A1 with case B, one can observe that it takes more time to fulfill the termination criteria if the heat losses are neglected. This is surprising because the temperature profiles at t2 and t3 and the vapor load at the top of the column are very similar. Again, the reason for the difference is the larger internal flow that leads to a faster dilution of impurities in the middle vessel. Process analysis was presented for a ternary mixture containing the same amounts of light-, intermediate-, and heavy-boiling components. We also carried out more simulations to investigate the influence of the feed composition using case C1 as the base set, which is typical for industrial-scale columns. The results show that purification of the middle-vessel product is the most time-consuming step for a wide variety of compositions. Conclusion Consideration of the entire startup offers realistic insight into the behavior of different column types, revealing analogous patterns. On the basis of a detailed process analysis, we can draw the following general conclusions: (a) Because adjustment of the distillate product purity was fast in the cases investigated, we recommend filling the distillate receiver as quickly as possible. Impurities of the middle boiler in the uppermost vessel can be easily diluted before the other key components reach specifications. Note that, in the case of a large excess of light-boiling component, it might be advantageous to start reflux before the holdup approaches the steadystate value. This is due to the larger time constant of the distillate receiver compared to the middle vessel in the mentioned case.
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008 823
However, temperature control can be used to avoid the accumulation of the middle-boiling component in the distillate receiver. (b) Because of condensation of the vapor during startup, the low-boiling component can accumulate in the middle vessel. The relative amount might be small. However, the impact on the process duration is high, particularly if high product purities are demanded. With respect to this observation, Gruetzmann et al.27 suggested bypassing the intermediate vessel during startup. Simulations presented in this contribution followed this recommendation consistently. (c) Moreover, we advise that the valve to the middle vessel be opened if most of low boiler is at the top of the column. This could be easily realized by a temperature sensor located in the upper part of the column used as an indicator. Because only the light and middle boilers reaches the upper column section, the boiling temperature is directly related to the composition. (d) We further recommend filling the middle vessel while providing a minimum reflux to the lower column section. The effect is twofold. Impurities trapped in the middle vessel will be diluted immediately, and separation of the middle and high boilers in the lower column section can proceed. Note that, because the presented conclusions mainly depend on the time constants of the product vessels, the process characteristics are independent of the relative volatilities for any ideal zeotropic mixture. Summary For the first time, a structured and extensive process analysis has led to general recommendations addressing improved operation of a middle-vessel batch distillation column. This work was divided into two major parts: modeling and simulation. A mathematical model that can handle the physical phenomena occurring during startup was presented. This dynamic model served as the basis for efficient simulation studies. It was clearly illustrated that use of a simplified initialization cannot lead to the same conclusions. The subsequent process analysis included relevant parameters related to the startup. As a result, we were able to provide general heuristics to improve the process for the first time. This is very important in the face of industrial applications of the process because previous works on this topic avoided making wide-ranging suggestions. Because the analysis revealed a number of optimization potentials, future work on dynamic optimization studies including startup is in progress. Nomenclature a ) specific surface area (m2 m-3) A ) surface area (m2) C ) constant used in eq 9 C1, C2 ) constants used in eq 11 cp ) specific heat capacity (J kg-1 K-1) d ) diameter (m) f ) generic function F ) vapor load (Pa0.5) g ) gravitational constant (m s-2) h ) specific enthalpy (J mol-1) H ) enthalpy (J) HU ) holdup (mol or L) k ) overall heat-transfer coefficient (W m-2 K-1) K ) equilibrium constant K1, K2 ) constants used in eq 14 L ) molar liquid flow (mol s-1) L ) length (m)
m ) mass (kg) N ) number of theoretical stages p ) pressure (Pa) q1 ) overall heat-transfer ratio (W W-1) q2 ) heating rate (K s-1) Q ) heat duty (W) t ) time (s) T ) temperature (K) u ) velocity (m s-1) V ) molar vapor flow (mol s-1) w ) mass fraction (kg/kg) x ) molar liquid fraction (mol mol-1) y ) molar vapor fraction (mol mol-1) Superscripts liq ) liquid vap ) vapor vl ) vapor-liquid, refers to vapor pressure Subscripts acc ) accumulation, refers to heat transfer to the equipment ambient ) refers to ambient conditions B ) bottom cond ) condenser D ) distillate dyn ) dynamic dyn0 ) refers to the dynamic liquid holdup without countercurrent flow equipment ) column internals and column shell F ) feed liq ) liquid loss ) refers to heat transfer to the surroundings m ) mean value max ) maximum i ) component index internals ) column internals, i.e., packings, accumulator, distributor j ) tray index p ) particle reb ) reboiler shell ) column shell stat ) static vap ) vapor Greek Symbols ) void fraction (m3 m-3) F ) density (mol m-3 or kg m-3) ν ) kinematic viscosity (m2 s-1) η ) dynamic viscosity (Pa s) σ ) surface tension (N m-1) τ ) hydraulic time constant (s) ∆ ) difference Ψ ) friction factor Dimensionless Numbers Bo ) Bond number Re ) Reynolds number
824
Ind. Eng. Chem. Res., Vol. 47, No. 3, 2008
Acknowledgment We gratefully acknowledge financial support from the MaxBuchner-Forschungsstiftung (MBFSt 2553). Literature Cited (1) Davidyan, A. G.; Kiva, V. N.; Meski, G. A.; Morari, M. Batch Distillation in a Column with a Middle Vessel. Chem. Eng. Sci. 1994, 49, 3033. (2) Barolo, M.; Botteon, F. Simple method of obtaining pure products by batch distillation. AIChE J. 1997, 43 (10), 2601. (3) Gruetzmann, S.; Fieg, G.; Kapala, Th. Theoretical analysis and operating behaviour of a middle vessel batch distillation with cyclic operation. Chem. Eng. Process. 2006, 45 (1), 46. (4) Warter, M.; Demicoli, D.; Stichlmair, J. Operation of a Batch Distillation Column with a Middle Vessel: Experimental Results for the Separation of Zeotropic and Azeotropic Mixtures. Chem. Eng. Process. 2004, 43, 263. (5) Hasebe, S.; Abdul Aziz, B.; Hashimoto, I.; Watanabe, T. Optimal Design and Operation of Complex Batch Distillation Column. In IFAC Workshop on Interaction between Process Design and Process Control; Pergamon Press: London, 1992; p 177. (6) Skogestad, S.; Wittgens, B.; Litto, R.; Sørensen, E. Multivessel Batch Distillation. AIChE J. 1997, 43 (4), 971. (7) Furlonge, H. I.; Pantelides, C. C., Sørensen, E. Optimal Operation of Multivessel Batch Distillation Columns. AIChE J. 1999, 45 (4), 781. (8) Low, K. H.; Sørensen, E. Simultaneous Optimal Configuration, Design and Operation of Batch Distillation. AIChE J. 2005, 581 (6), 1700. (9) Wittgens, B.; Litto, R.; Sørensen, E.; Skogestad, S. Total Reflux Operation of Multivessel Batch Distillation. Comput. Chem. Eng. 1996, S20, 1041. (10) Gruetzmann, S.; Temprano, M.; Fieg, G.; Kapala, Th. Fu¨hrungskonzepte zur Minimierung der Prozessdauer der zyklischen BatchRektifikation. Chem. Ing. Technol. 2005, 77 (8), 1101. (11) Barolo, M.; Guarise, G.B.; Rienzi, A.; Trotta, A., Nonlinear ModelBased Startup and Operation Control of a Distillation Column: An Experimental Study. Ind. Eng. Chem. Res. 1994, 33, 3160. (12) Fieg, G.; Wozny, G.; Kruse, C. Experimental and Theoretical Studies of the Dynamics of Startup and Product Switchover Operations of Distillation Columns. Chem. Eng. Process. 1993, 32, 283. (13) Ganguly, S.; Saraf, D. N. Startup of a Distillation Column Using Nonlinear Analytical Model Predictive Control. Ind. Eng. Chem. Res. 1993, 32, 1667. (14) Ruiz, C. A.; Cameron, I. T.; Gani, R. A Generalized Dynamic Model for Distillation Columns III. Study of Startup Operations. Comput. Chem. Eng. 1988, 12, 1.
(15) Wang, L.; Li, P.; Wozny, G.; Wang, S. A Startup Model for Simulation of Batch Distillation Starting from a Cold State. Comput. Chem. Eng. 2003, 27, 1485. (16) Elgue, S.; Prat, L.; Cabassud, M.; Le Lann, J. M.; C’ezerac, J. Dynamic Models for Start-up Operations of Batch Distillation Columns with Experimental Validation. Comput. Chem. Eng. 2004, 28, 2735. (17) Gruetzmann, S.; Fieg, G.; Kapala, Th. Dynamic Modeling of Complex Batch Distillation Starting from Ambient Conditions. In 16th European Symposium on Computer Aided Process Engineering and 9th International Symposium on Process Systems Engineering; Marquardt, W., Pantelides, C., Eds.; Elsevier: Amsterdam, 2006; pp 865-870. (18) Kruse, Ch. Theoretische und experimentelle Entwicklung zeitoptimaler Anfahr- und Produktwechselstrategien fu¨r die Rektifikation. Fortschr.Ber. VDI, Reihe 3 1995, (420). (19) Flender, M. Zeitoptimale Strategien fu¨r Anfahr- und Produktwechselvorga¨nge an Rektifizieranlagen. Fortschr.-Ber. VDI, Reihe 3 1999, (610). (20) Engel, 2004.; Stichlmair, J.; Geipel, W. Fluid Dynamics of Packings for Gas-Liquid Contactors. Chem. Eng. Technol. 2001, 24 (5), 549. (21) Mackowiak, J. Fluiddynamik Von Fu¨llko¨rpern und Packungen; Springer-Verlag: Berlin, 2003. (22) Ascher, U. M.; Petzold, L. Computer Methods for Ordinary Differential Equations and Differential Algebraic Equations; SIAM: Philadelphia, PA, 1998. (23) Pantelides, C. C.; Gritsis, D.; Morison, K. R.; Sargent, R. W. H. The Mathematical Modelling of Transient Systems Using DifferentialAlgebraic Equations. Comput. Chem. Eng. 1988, 12 (5), 449. (24) Fa´bia´n, G.; van Beek, D. A.; Rooda, J. E. Index Reduction and Discontinuity Handling Using Substitute Equations. Math. Comput. Modell. Dyn. Syst. 2001, 7 (2), 173. (25) Rix, A. Modellierung und Prozessfu¨hrung wa¨rmeintegrierter Destillationskolonnen. Fortschr.-Ber. VDI, Reihe 3 1995, (549). (26) Aspen Custom Modeler Release 2004.1, Aspen Modeler Reference Guide; Aspen Technology, Inc.: Cambridge, MA, 2005. (27) Gruetzmann, S.; Fieg, G. Improvements in Configuration of a Middle Vessel Batch Distillation by Means of Rigorous Dynamic Modelling and Simulation. In Proceedings of the 15th IASTED International Conference on Applied Simulation and Modelling; Hamza, M. H., Ed.; ACTA Press: Calgary, Canada, 2006; pp 106-110.
ReceiVed for reView May 11, 2007 ReVised manuscript receiVed July 25, 2007 Accepted October 30, 2007 IE070667V