**THERMODYNAMICS AND SEPARATION PROCESSES**

**Dynamic simulation of flash drums using rigorous physical property calculations**

**F. M. Gonçalves; M. Castier ^{*}; O. Q. F. Araújo**

Escola de Química, Universidade Federal do Rio de Janeiro, Fax: +(55) (21) 2562-7567, Cx P. 68542, CEP: 21949-900, Rio de Janeiro - RJ, Brazil. E-mail: castier@eq.ufrj.br

**ABSTRACT**

__N__flash). A new aspect of our dynamic simulations is the use of direct iterations in phase volumes (instead of pressure) for solving the algebraic equations. It was also found that an iterative procedure previously suggested in the literature for UV

__N__flashes becomes unreliable close to phase boundaries and a new alternative is proposed. Another unusual aspect of this work is that the model expressions, including the physical properties and their analytical derivatives, were quickly implemented using computer algebra.

**Keywords:** Dynamic simulation; Flash; Equation of state; Phase equilibrium; Hydrocarbons; Differential-algebraic equations.

**INTRODUCTION**

Vessels can be the most hazardous pieces of equipment in a chemical plant and, therefore, knowledge of their dynamic behavior is important for proper design and operation (Driedger, 2000). Nevertheless, typical approaches for their dynamic simulation imbed simplifications such as ideal liquid and/or vapor phases and temperature-independent heat capacities, inappropriate for describing applications at high pressures that may require the use of equations of state (EOS) for more rigorous calculations. Recent formulations are shown in Table 1, in which some features of this work are included for comparison. In most cases, a flash with specified values of internal energy, volume and mole numbers (UV__N__ flash) is at the core of the calculations executed at each time step. Müller and Marquardt (1997) and Gopal and Biegler (1997) modeled the liquid phases in their examples using excess Gibbs free energy (G^{E}) models, which, from the computational point of view, are much simpler than EOS. Furthermore, Gopal and Biegler (1997) assume perfect pressure control and specify the temperature trajectory in the drum, what results in a much simpler version of the dynamic flash problem. The articles of Saha and Carroll (1997) and Michelsen (1999) are based on UV__N__ flash formulations that use nested-loops since pressure is an iteration variable. This has the disadvantage of requiring the solution of the EOS to obtain the molar volume of each phase in every iteration. An important suggestion made by Michelsen (1999), but not pursued in his work, was to investigate UV__N__ flashes with direct iterations in temperature, phase volumes and moles numbers of each component in each phase. This is the approach used herein.

In this work, we simulate the dynamic behavior of storage tanks and flash drums. A new aspect of these dynamic simulations is the use of direct iterations in phase volumes (instead of pressure) for solving the algebraic equations. It was also found that an iterative procedure previously suggested in the literature for UV__N__ flashes (Michelsen, 1999) becomes unreliable close to saturation points (e.g., bubble or dew points) and a new alternative is proposed. Another unusual aspect of this work is that all the model expressions, including the thermodynamic properties and their analytical derivatives, were quickly implemented using computer algebra.

**FORMULATION**

We consider a drum with input streams and S_{W} output streams. The drum volume (V) is constant and known. The energy and mass balance equations are:

In these equations, denotes time, U is the internal energy, and N_{i} is the number of moles of component i in the tank. The molar flowrate of component i and the enthalpy flowrate of input stream j are denoted by and , respectively. The symbols and denote analogous values for output stream. The heat load in the drum is denoted by and n_{c} is the number of components. Integration of the (n_{c} + 1) differential equations (1) and (2) gives the evolution of U and __N__. Therefore, U, V, and __N__ are known at any time step. Even though the algebraic and differential equations can be uncoupled in some cases, for generality, they are numerically integrated in our procedure.

We assume thermodynamic equilibrium in the drum at any time and that the intensive properties of the output streams are equal to those of the phase in the drum from which they are withdrawn. For specified UV__N__ values, equilibrium can be calculated by maximizing the entropy, but this is cumbersome since most EOS in practical use have the form P=P(T,V,__N__). We followed Michelsen (1999), using a function Q_{F} given by:

where A is the Helmholtz free energy, U* is the specified internal energy (the value of U at a given time), and R is the universal gas constant. In the case of vapor-liquid flash problems, Michelsen suggested the differentiation of Q_{F} with respect to the number of moles of each component in one of the phases and to the logarithms of temperature and of one of the phase volumes (e.g., the vapor phase), Using these derivatives, we observed that convergence of the UV__N__ flash was very difficult or could not be achieved when the amount of one of the phases was small (i.e., close to saturation points, such as bubble or dew points). One of the new findings of this work is the elimination of this difficulty by using derivatives with respect to temperature and phase volumes. Therefore, for a system containing phases, we used the following set of algebraic equations constituted by derivatives of Q_{F}:

where U_{k} and P_{k} are the internal energy and pressure in phase; µ_{ij} is the chemical potential of component in phase, and denotes the phase in which the largest amount of component is located. The problem was formulated assuming that the volume of one of the phases (phase K) is a dependent variable, calculated from the volume conservation equation. For numerical convenience, the largest phase volume at each iteration was assumed to be a dependent variable. Similarly, the largest number of moles in a phase for each component was assumed to be a dependent variable.

Equations 4-6 represent a set of (n_{p}+n_{c}(n_{p}-1))nonlinear algebraic equations. In the case of a single-phase system, the set reduces to Eq. (4). One of the advantages of this formulation is that the Jacobian matrix of the set of equations is the Hessian matrix ofQ_{F}; therefore, this matrix is symmetrical and only of its terms need to be evaluated. The full Newton-Raphson step at each iteration is obtained by solving the set of linear equations:

where __D__V and __D__N are vectors containing the independent volumes and component mole numbers arranged in a convenient order. The full Newton-Raphson step occasionally leads to non-physical values such as negative volumes or mole numbers. In this case, the step size was reduced but the search direction was preserved. The analytical derivatives needed to compute the Jacobian matrix are presented by Gonçalves (2001).

*ode15s*available in MATLAB was initially used for solving the resulting non-linear DAE system. In our preliminary tests, however, this procedure failed and reported DAE index larger than one, a point we are investigating further. Our approach was then to use the solver ability to handle stiff systems while declining to use its DAE solving capability. In the second implementation, the differential equations of the model (Eqs. (1) and (2)) were integrated using the Bulirsch-Stoer algorithm as implemented by Press et al. (1992). The initial step size was set equal to 1/1000 of the total time span. The maximum rate of increase in time step between consecutive steps was set equal to 1/0.999. It should be remarked that the formulation is general and other integration methods can be used.

For the calculation of the physical properties of the input streams, which may vary in time and contain one or more phases, a TP__N__ flash procedure (Michelsen, 1982b) was used. The solution of a TV__N__ flash (Espósito et al., 2000) provides the initial drum condition. The algebraic UV__N__ flash equations usually converged in one to three iterations, except close to saturation points, where four to eight iterations were typically required. A given phase j was removed whenever was smaller than 10^{-6}. In this case, the amount of the removed phase was added to the other phase and the algebraic equations were changed accordingly.

For determining whether a new phase should be added to the system, the global stability (Michelsen, 1982a) test was used. The mole fractions and molar volume of the incipient phase in the global minimum of the stability function were used as initial estimates for the composition of the new phase in the drum. The new phase was initialized with % of the total number of moles present in the drum.

The solution of the UV__N__ flash corresponds to a saddle point of the Q_{F}-function. Similarly to previous results reported in the literature (Michelsen, 1987), we had difficulties to converge three-phase problems using a saddle-point formulation. We are currently investigating the numerical behavior of the UV__N__ flash formulated as a minimization problem, by adding a penalty term to the Q_{F}-function that depends on the squared deviation between the calculated and specified values of the internal energy.

Physical properties were computed using the Peng-Robinson (1976) EOS with the one-fluid van der Waals mixing rules. Pure component properties were taken from Reid et al. (1987). Binary interaction parameters were set equal to 0 (zero). Several properties are necessary for solving the UV__N__ flash: internal energy, pressure, chemical potential (or the logarithm of the fugacity) of each component, and their derivatives with respect to temperature, volume and mole numbers. It should be noted that derivatives at constant phase volumes are not usually available, in direct form, in standard physical property packages. All the necessary expressions were quickly obtained and implemented using the *Thermath* computer algebra package (Castier, 1999). Even though these physical properties and their derivatives all have real values, intermediate variables may have complex values, as consequence of the direct iterations on temperature, volumes (instead of pressure), and mole numbers, and because the intermediate variables are automatically defined by the computer algebra package. For this reason, complex arithmetic was used to calculate these intermediate variables.

**RESULTS**

The examples illustrate several types of specifications, such as pure components or mixtures, and closed or open systems (only with input or output streams, or with both). In some of the examples, the number of phases (one or two) remains constant throughout the simulation. In others, there is appearance or disappearance of a phase, or both.

**Example 1**

In this example, a new phase appears during the simulation. We consider a tank with a volume of 30 m^{3} being filled with pure methane. We solve two problems, 1a and 1b, whose specifications are in Table 2. In both cases, the specified heat load corresponds to a heat removal whose rate is time-dependent. Figures 1 and 2 show the evolution of temperature and fluid pressure in the tank in problem 1a. With this set of specifications, there is no condensation of methane in the tank. Figures 3 and 4 present the evolution of the same two variables in problem 1b, showing that both variables pass through a maximum and then decrease after the beginning of heat removal from the tank. At approximately 829 min, there is the onset of a second phase in the tank, what appears in Figs. 3 and 4 as sudden changes in slope.

]]>

**Example 2**

In this example, a phase disappears during the simulation. We consider a mixture of ethane(1) + propene(2) + propane(3) + isobutane(4) + n-butane(5) + n-pentane(6) whose composition is in the range of liquefied petroleum gases (LPG). At the initial condition, presented in Table 3, there are a liquid and a vapor phase in the drum. A vapor stream is continuously withdrawn from the drum and a constant heat load is provided to prevent a large temperature decrease. The problem specification is in Table 4. Figures 5 and 6 show the evolution of temperature and pressure in the drum. Figure 6 is especially interesting because it shows that the pressure goes through a maximum in the two-phase region at approximately 220 min. The moment when the liquid phase disappears from the drum, at approximately 439 min, is a local minimum in pressure with a small cusp, and there is a second maximum in pressure, at approximately 500 min.

]]>

]]>

**Example 3**

We use the same mixture and initial conditions of Example 2, shown in Table 3. In this case, however, there is one input stream and two output streams, corresponding to the withdrawal of the liquid and the vapor phases. The problem specification is presented in Table 5. Given that the flowrate of the output streams is smaller than that of the input stream, material accumulates in the drum. Examining the output specifications, the flowrate of the liquid output is smaller than that of the vapor output. For this reason, the heavier components accumulate at a faster rate, as can be observed in Figure 7. The temperature (Fig. 8) passes through a maximum very close to the beginning of the simulation time and then through a minimum at approximately 70 min. The pressure (Fig. 9) passes through a minimum at approximately 24 min. In this problem, there is no need for phase addition or removal because there are two phases in drum throughout the simulation time.

]]>

**Example 4**

In this example, we again use the mixture and initial conditions of Example 2, presented in Table 3. In this example, the tank is closed, but there is a heat load that varies with time according to:

(T in min and Q in kJ/min)

The system was simulated between 0 min and 400p=1256.64 min, corresponding to two complete disturbance cycles. At the initial condition, there are a liquid and a vapor phase in the drum. Figure 10 shows the evolution of the vapor phase volume. Where a horizontal line represents this volume, there is only the vapor phase in drum. In Figure 11, we show the behavior of the mole fraction of two of the components, propane and n-pentane, in the liquid phase. Their mole fraction evolutions have opposite patterns because propane is one the lighter components while n-pentane is the heaviest component in the mixture. The interruptions in the lines correspond to the time intervals when there is no liquid phase in the drum. In Figure 12, we show the temperature evolution, in which the discontinuities in the derivative of the temperature with respect to time occur at the points of phase appearance or disappearance. The pressure in the system oscillates between 0.4633 and 0.5155 MPa. This example shows a situation in which there are both phase appearance and disappearance in a single simulation. Given the specified disturbance, which is periodic, it should be expected that the system returned to the initial condition after each complete cycle, which lasts 200p min. This was indeed observed for all variables, showing the consistency of the calculated results.

]]>

]]>

**CONCLUSIONS**

The dynamic behavior of flash drums was simulated using a formulation adequate for phase modeling with equations of state (EOS). A new aspect of our dynamic simulations is the use of direct iterations in phase volumes (instead of pressure) for solving the algebraic equations, what has the advantage of avoiding the need for nested loops. It was also observed that an iterative procedure for the UV__N__ flash problem, suggested in the literature, based on the use of logarithms of temperature and phase volumes becomes unreliable close to phase boundaries. We advocate the use of temperature and phase volumes as iterated variables because of the better convergence behavior.

**ACKNOWLEDGEMENTS**

F.M.G. acknowledges a scholarship from Agência Nacional do Petróleo (ANP, Brazil). The authors are grateful to ANP, CNPq/CTPetro/Brazil, FAPERJ and PRONEX (Contract 124/96) for their financial support.

**NOMENCLATURE**

**Latin Letters**

A | Helmholtz free energy | (-) |

f_{ij} | molar flowrate of component in input stream | j |

enthalpy flowrate of input stream | j | |

enthalpy flowrate of output stream | k | |

K_{i} | distribution coefficient of component | i |

n_{c} | number of components | (-) |

N_{ij} | number of moles of component in phase | j |

N_{i} | number of moles of component in the drum | (-) |

n_{P} | number of phases | (-) |

P | pressure | (-) |

heat load | (-) | |

Q_{F} | flash function | (-) |

R | universal gas constant | (-) |

S_{f} | number of input streams | (-) |

S_{W} | number of output streams | (-) |

t | time | (-) |

T | temperature | (-) |

U | internal energy | (-) |

V | drum volume | (-) |

W_{ik} | molar flowrate of component in output stream | k |

**Greek Letter**

µ_{ij} | chemical potential of component in phase | j |

**REFERENCES**

Castier, M., Automatic Implementation of Thermodynamic Models Using Computer Algebra, Computers and Chemical Engineering, 23, 1229 (1999). [ Links ]

Driedger, W.C., Controlling Vessels and Tanks, Hydrocarbon Processing, 79, No. 3, 101 (2000). [ Links ]

Espósito, R.O., Castier, M., and Tavares, F.W., Calculations of Thermodynamic Equilibrium in Systems Subject to External Fields, Chem. Engng Sci., 55, 3495 (2000). [ Links ]

Gonçalves, F.M., M.Sc. Thesis, Federal University of Rio de Janeiro, Brazil (2001). [ Links ]

Gopal, V. and Biegler, L.T., Nonsmooth Dynamic Simulation with Linear Programming Based Methods, Computers and Chemical Engineering, 21, 675 (1997). [ Links ]

Michelsen, M.L., The Isothermal Flash Problem. I. Stability Analysis, Fluid Phase Equilibria, 8, 1 (1982a). [ Links ]

Michelsen, M.L., The Isothermal Flash Problem. II. Phase Split Calculations, Fluid Phase Equilibria, 8, 21 (1982b). [ Links ]

Michelsen, M.L., Multiphase Isenthalpic and Isentropic Flash Algorithms, Fluid Phase Equilibria, 33, 13 (1987). [ Links ]

Michelsen, M.L., State Function Based Flash Specifications, Fluid Phase Equilibria, 158-160, 617 (1999). [ Links ]

Müller, D. and Marquardt, W., Dynamic Multiple-Phase Flash Simulation: Global Stability Analysis versus Quick Phase Determination, Computers and Chemical Engineering, 97, S817 (1997). [ Links ]

Peng, D.Y. and Robinson, D.B., A New Two-Constant Equation of State, Ind. Eng. Chem. Fundamentals, 15, 59 (1976). [ Links ]

Press, W.H., Teukolsky, S.A. Vetterling, W.T., and Flannery, B.P., Numerical Recipes in FORTRAN, 2^{nd} ed., Cambridge University Press, Cambridge (1992). [ Links ]

Reid, R.C., Prausnitz, J.M., and Poling, B.E., The Properties of Gases and Liquids, 4^{th} ed., McGraw-Hill, New York (1987). [ Links ]

Saha, S. and Carroll, J.J., The Isoenergetic-Isochoric Flash, Fluid Phase Equilibria 138, 23 (1997). [ Links ]

(Received: October 25, 2003 ; Accepted: November 10, 2005)

]]>

* To whom correspondence should be addressed

]]>