Differential-Algebraic Approach to Dynamic Simulations of Flash Drums with Rigorous Evaluation of Physical Properties

Résumé — Simulation dynamique de séparateurs de mélanges liquide-vapeur avec une approche différentielle algébrique et évaluation rigoureuse des propriétés physiques — La dynamique de séparateurs de mélanges liquide-vapeur est simulé avec des propriétés physiques rigoureusement calculées en utilisant une équation d’état pour modéliser chaque phase. La formulation fournit un ensemble d’équations différentielles algébriques (DAE), dont les équations différentielles décrivent les bilans matérielles et énergétiques et les équations algébriques décrivent les conditions d’équilibre à l’intérieur du séparateur. PSIDE (Parallel Software for Implicit Differential Equations) [Lioen et al. , 1998] est utilisé pour résoudre l’ensemble des DAE avec efficacité. Dans cette approche, les équations sont résolues simultanément, avec des itérations directes de température, de volumes de phase et de nombre de moles de chaque composant dans chaque phase. Les résultats montrent l’efficacité de cette méthode de simulation des séparateurs de mélanges liquide-vapeur. Abstract — Differential-Algebraic Approach to Dynamic Simulations of Flash Drums with Rigorous Evaluation of Physical Properties — The dynamics of flash drums is simulated with rigorous physical properties calculations using


INTRODUCTION
The ability to predict the dynamic behavior of separation drums is important to devise control strategies that ensure product quality and safe operation, especially at high pressures. Quite often, dynamic simulations of flash drums are based on simple thermodynamic models, such as ideal vapor and liquid phases. To account for deviations from ideal liquid phase behavior, some articles [1,2] include activity coefficients in their models, but even with this extension, the pressure range in which the model can be used with confidence is limited. Calculations of heat effects are frequently based on temperature-independent heat capacities, what is certainly a good assumption when the temperature variation during the simulated time interval is not very large. Though, for larger temperature variations, use of temperature-dependent heat capacities will typically improve accuracy. Furthermore, in some cases, the time evolution of pressure is assumed beforehand, based on the assumption it is controlled. More realistic simulations, however, should be able to predict how pressure evolves. Modern equations of state (EOS) can help overcome many of these limitations as they can model the phase behavior of many systems, including mixtures containing polar components, polymers, and electrolytes, up to high pressures. Therefore, it is interesting to develop a framework that can be used to simulate the dynamics of separation vessels with phases modeled by EOS. Few papers on vessel dynamics discuss the use of EOS [3][4][5]. This is by no means a trivial extension of dynamic models based on simpler assumptions about phase behavior because: -vapor-liquid equilibrium calculations tend to be more difficult at high pressures than at low pressures, because of the larger departure from ideal phase behavior and the proximity of the mixture critical point; -depending on calculation algorithm and problem specifications, it may be necessary to solve the equation of state to find the molar volume at a given condition of temperature, pressure and system composition. This is, by itself, an iterative problem embedded inside the general structure of the simulation model; -expressions for properties, such as fugacity coefficients, derived from EOS commonly used for process simulation usually have more terms than those derived from excess Gibbs free energy models, such as activity coefficients. Furthermore, in the EOS most commonly used for chemical process design, properties are explicit functions of temperature, volume, and number of moles of each component in a given mixture. Depending on the values used for these variables, EOS may predict, for example, negative pressures and negative fugacity coefficients, both of which have no physical meaning. The model developed here assumes phase equilibrium inside the vessel at each instant throughout the simulation. With this assumption, the model constitutes a system of differential algebraic equations (DAE). The mass and energy balances for the separation vessel are differential equations that enable computations of the time evolution of internal energy (U) and number of moles of each component of the fluid inside the separator (N). Since the tank volume (V) is constant, the values of U and N should be determined through numerical integration of mass and energy balances at each instant. The algebraic equations of the model come from the phase equilibrium conditions at specified values of UVN, the so-called UVN flash problem, which is at the core of our simulations of flash drum dynamics.
In spite of its importance for separator dynamics, few papers address the UVN flash problem. The conditions for thermodynamic equilibrium, at specified UVN, maximize the system's entropy. The difficulty with the direct application of this formulation is the need for partial derivatives of entropy with respect to internal energy, volume, and number of moles. However, the most commonly thermodynamic models used for chemical process design are written as functions of temperature, volume, and number of moles. A family of functions has been proposed [4] to overcome this difficulty for many specifications, and in particular, used a UVN flash algorithm with "nested" loops, in which pressure is an iteration variable. A disadvantage of this approach is the need to solve the EOS for each phase in every iteration, since most EOS have the form P = P(T, V, N) . The direct use of temperature, phase volumes, and mole number of each component in each phase as iteration variables has been suggested [4]. This approach was adopted [6] in the development of a procedure to simulate flash drum dynamics, in which a conventional integrator for ordinary differential equations was used and the algebraic equations of the model were solved separately in every step. Attempts to solve the system directly using a procedure available in MATLAB for DAE (ode15s solver) were unsuccessful [6]. The most likely reason is that the ode15s integrator can only solve linearly implicit ODE systems of general form: M(t) · y' = f(t, y) with a mass matrix M(t) that is non-singular, (usually) sparse and exclusively time dependent [7]. The proposed formulation led to index-1 and index-2 DAE systems that can be rearranged in that form only after an unacceptable manipulation of the algebraic constraints. Lee et al. [5] simulated the dynamic behavior of flash drums using the Soave-Redlich-Kwong equation of state to calculate physical properties. This DAE problem was solved using MATLAB with application of the recursive equation error method to overcome the difficulties introduced by the use of the EOS.
In this work, we take a different path that enables us to solve the DAE system directly and simultaneously. A significant decrease on computational cost has been verified using this procedure in comparison with previous work [6], in which all the algebraic equations were solved separately in every step. Using PSIDE, the algebraic constraints were solved together with the discretized form of the differential equations by a very efficient Newton process. Additionally, PSIDE is capable of solving DAE systems with differential index up to 3. It is important to notice that typical software used to solve ODE and DAE systems, such as DASSL [8], are able to solve systems of index up to 1 only.
Hence, the main purpose of this work is to predict the dynamics of flash drums, with rigorous calculations of physical properties using EOS and solving the set of DAE's directly. The next section presents the PSIDE features. We then discuss the model and present results of its application. The article ends with the conclusions drawn from this research.

PSIDE
We used PSIDE [9] to integrate the DAE system. PSIDE is a code written in FORTRAN-77 capable of solving implicit differential equations of the form: f(t, y, y') = 0, f, y ∈ ℜ n , t 0 ≤ t ≤ t end , y(t 0 ) = y 0 and f(t 0 , y 0 , y' 0 ) = 0. PSIDE integrator is the four-stage Radau IIA method. The nonlinear systems are solved by a modified Newton process, in which every Newton iterate itself is computed by means of a parallel iterative linear system solver for Runge-Kutta [10]. The Radau IIA method is a class of implicit Runge-Kutta method (IRK), powerful for the numerical solution of implicit differential equations. As a Radau IIA method, PSIDE combines high order with the property of L-stability. It is a robust code that allows solving algebraic-differential systems with differential index up to 3, unlike classical software that generally have problems when solving DAE systems with index higher than 1.
The differential index is, in a general way, closely linked to the number of times that the system or a part of it should be analytically differentiated in order to obtain a purely differential system through algebraic manipulations. Therefore, the procedure of differential index reduction is quite costly and the possibility of working directly with higher-index systems represents a great advantage.
The code allows the user to supply the analytical Jacobian of the system, but, if not provided, it is calculated numerically, using finite differences. The user must provide the socalled mass matrix (M), given by: where g (t, t, y') are the equations of the system and y' represents the derivatives of the dependent variables y with respect to the independent variable t.

FORMULATION
Consider the flash drum of constant volume shown in Figure 1. The material balance of each component is given by: where t denotes time, N i is the mole number of component i in the drum, z i is the mole fraction in the input stream, x i and y i are the mole fractions of component i in the liquid and in the vapor outputs, respectively, and c is the number of components in the system. In Equation 1, F, L, and G denote the molar flowrates of the feed stream, liquid output, and vapor output, respectively. Equation 1 can be readily adapted to problems with fewer or more input and output streams. Example 1, in the next section, will illustrate a situation with single input and output streams and a single phase inside the drum.
Assuming that the mole fractions of the output streams are identical to those of the phases inside the drum from where they are withdrawn, then: where N i,1 and N i,2 are the number of moles of component i in the liquid and vapor phases inside the drum, respectively.
. . The global energy balance is: where H · F , H · L and H · G are the enthalpy flowrates of the input, liquid output, and vapor output streams, respectively; Q · is the heat load in the drum, and U is the internal energy of the system. To calculate the enthalpies of each phase or stream, we compute the ideal gas contribution using heat capacity at constant pressure data for ideal gases (Reid et al. [11]), to which we add the residual enthalpy, calculated according to the EOS. The internal energies of the phases inside the drum are calculated in similar fashion by using ideal gas heat capacity at constant volume data and residual internal energy from the EOS. Equations 2 and 3 constitute a set of (c + 1) ordinary differential equations (ODE) that enables the computation of N and U as functions of time. At specified values of UVN, the equilibrium condition maximizes the system's entropy. To overcome the difficulty that the internal energy is not an independent variable in the thermodynamic models most commonly used for process design, Michelsen [4] defined a new function, Q F , given by: (4) where A denotes the Helmholtz free energy, U is the internal energy in each step (obtained by integration of the energy balance, Equation 3), R is the universal gas constant, and T is the temperature. The equilibrium condition occurs at a stationary point of Q F , but unfortunately not a minimum since = -Cp/T < 0, and the Hessian matrix is therefore not positive definite, as required for a minimum. In this case, the stationary point is a saddle point, where its value is equal to the system's entropy maximum. The advantage of this formulation is that the independent variables of function Q F are T, V, and N that are also the independent variables of the EOS most commonly used for chemical process design. The Jacobian matrix of this set of equations is the Hessian matrix (matrix of which the entries are the second partial derivatives of the function) of Q F . An important feature of the Hessian matrix is that it is symmetrical, reducing by approximately half the effort for calculating derivatives.
One should notice that, because of volume conservation and mass balances, not all phase volumes and mole numbers can vary independently: the value in one of the phases can be taken as a dependent variable. Differentiation of function Q F with respect to the independent variables in a system with p phases gives [6]: i = 1, ..., c j = 1, ..., p j ≠ J In these equations, U j and P j are the internal energy and the pressure of the phase j, respectively, μ ij and f ij are the chemical potential and fugacity of component i in phase j, respectively. Subscript J denotes the dependent phase. Equation 5 ensures that the system internal energy matches the value specified at each time step. Equations 6 and 7 provide the conditions for equal phase pressures and chemical potentials. Volume conservation and material balances provide additional equations: where V j is the volume of phase j. namely: N i , N ij , V j , P j , T and U, with the component and phase subscripts running from 1 to c and 1 to p, respectively.
To characterize the differential index of the resulting DAE system we can identify 3 kinds of variables: -Differential variables, i.e., variables submitted to time derivatives: N i and U; -Algebraic variables, i.e., variables not submitted to time derivatives: P j , V j , T, L, and G. Algebraic constraints contain variables, P j , V j , and T. Although, output flowrates L and G do not appear in any algebraic equation. The first kind of algebraic variables (P j , V j and T) were classified as explicit algebraic variables and the second kind (L and G) as non-explicit algebraic variables. Based on this variable classification, we can present a general model that contains all the properties of separation drums dynamic mathematical model. The proposed model is described by the equations: , constrained by the algebraic equations: g[t, x(t), y(t)] = 0. To obtain the time derivative of variable y, this algebraic constraint must be differentiated, resulting in: , where g y ≠ 0. This procedure shows that the system is at least an index-1 DAE system. As no information about the algebraic variable z is provided, a new equation must be specified in the problem.

Two possibilities are considered: a) A new algebraic constraint h[t, x(t), y(t), z(t)]
= 0 is specified, where variable z(t) appears. To obtain the time derivative of variable z, this algebraic constraint must be differentiated, resulting in: In this case, the index of the DAE system remains equal to 1, and the algebraic variables y and z can be both classified as explicit algebraic variables.

b) A new algebraic constraint h[t, x(t), y(t)] = 0 is specified,
where variable z(t) still does not appear. In this case, to obtain the time derivative of variable z, this algebraic constraint must be differentiated twice. The first time differentiation provides: x, y, z] = 0, another time differentiation results in: In this case the system differential index is 2, algebraic variable y is an explicit algebraic variable and z is a nonexplicit algebraic variable. In Examples 1 to 4 below, we assume that the input flowrate, mole fractions, temperature and pressure are known, and that the output flowrates L and G are controlled, both being explicit algebraic variables (similar to case (a) above). In Example 5, we consider that the pressure is known (explicit algebraic variable), the output liquid flowrate L is controlled (explicit algebraic variable) and the output vapor flowrate G is unknown (non-explicit algebraic variable). In Example 6, we consider that mole fraction of one of the components in vapor phase is specified (related with a differential variable), the output liquid flowrate L is controlled (explicit algebraic variable) and the output vapor flowrate G is unknown (non-explicit algebraic variable). DASSL and PSIDE codes have been applied to solve Examples 1 to 4. These examples are all index-1 DAE system. Examples 5 and 6 are both index-2 DAE system (similar to case (b) above) and only PSIDE code was capable of solving them, without using any reduction index procedure.
A TPN flash is solved to establish the thermodynamic state of the input stream(s) and the corresponding enthalpy flowrate(s) (H f ). Solution of a TVN flash [12] provided selfconsistent initial drum conditions. Thermodynamic properties were calculated using the Peng-Robinson [13] EOS with one-fluid van der Waals mixing rules. All the necessary thermodynamic expressions were obtained and implemented using the Thermath computer algebra package [14]. Pure component properties were taken from reference [11] and binary interaction parameters were set equal to 0 (zero).
The resulting model was implemented in FORTRAN and directly solved as a DAE system using the PSIDE routine.

RESULTS
The following examples illustrate the proposed approach for different specifications. In all of them, CPU time was below one minute in a 3 GHz Pentium IV processor, with 1 Gb of RAM, using an integration step Δt = 0.1 units of time.

Example 1
This first example [3] consists in filling a vessel with pure nitrogen. The vessel has inlet and outlet connections. There is a control valve in the outlet, which is initially closed. The coefficient of this valve is chosen to be equal to 2, the internal diameter of the pipe is 20 cm and the nominal valve size is 15 cm. This control valve is modeled according to reference [15]. It is also assumed that there is a linear heat load between the fluid and the wall of the vessel: Q · = -3 (t -30) [kJ/s]. Table 1 presents the specifications and initial conditions for this example. Figure 2 shows the dynamic change in the molar volume of N 2 and Figure 3 presents the evolution of the vessel temperature and pressure. As expected, the pressure increases with time. The initial heat load (for t < 30 s) is positive and becomes negative afterwards, causing a maximum in temperature.

Example 2
This example deals with a mixture of ethane(1), propene(2), propane(3), n-butane(4), isobutane(5) and n-pentane(6) whose composition is in the range of liquefied petroleum gases (LPG) [6]. Table 2 presents the initial conditions, in which a liquid and a vapor phase are in equilibrium inside the drum. There is a single input stream and a liquid and a vapor output streams. Table 3 presents the problem specification.
Mass accumulates in the drum because the total flowrate of the output streams is smaller than that of the input stream. Moreover, the heavier components accumulate at a faster rate because the liquid output flowrate is smaller than the vapor output flowrate, as shown in Figure 4. Figures 5a and 5b show mole fractions in the liquid and vapor phases, respectively, as functions of time. Figure 6 shows the evolution of temperature and pressure in the drum. The temperature profile is particularly interesting because it passes through a maximum and a minimum point as result of combining several effects. The input stream is hotter (300 K) than the mixture initially inside the drum (298.15 K). This has the immediate effect of increasing the drum temperature in the beginning of the simulation. This is later compensated by the effect of input stream expansion when it enters the tank, which is at lower pressure, causing a decrease in temperature.   The difference between the input stream and tank pressures becomes smaller afterward and the temperature increases as result of mass accumulation in the tank. Our results have a perfect match with those of the literature [6].

Example 3
We use the same mixture and initial conditions of Example 2, presented in Table 2. However, the tank is closed and there is a heat load that varies with time according to Q · 6.26sin (0.01t) [t in min and Q · in kJ/min].
The system was simulated during a period corresponding to two complete disturbance cycles (400π = 1256.64 min). At the initial condition, there are a liquid and a vapor phase in the drum. Figure 7 shows the evolution of the liquid phase volume in the tank. The heat load was chosen in such a way that the liquid phase almost disappears in each cycle. The liquid phase mole fractions of propane and n-pentane are shown in Figure 8 and they also present periodic behavior, with opposite patterns. As n-pentane is the heaviest component in the system, its mole fraction in the liquid phase tends to increase when the temperature increases. Propane is much lighter and its liquid phase mole fraction has opposite pattern.
The specified disturbance is periodic. Hence it should be expected that the system returned to the initial condition after each complete cycle. This was indeed observed for all variables, showing the consistency of the results presented here.

Example 4
In this example, we consider the same system as in Examples 2 and 3 and the tank is closed, as in Example 3. Here, to test  the consistency of our results, we compare a sinusoidal heat load Q · = sin (0.01t) with a piecewise heat load, as described in Table 4. We observe that the heat load defined in the fourth time range of Table 4 is the negative of that one defined in the second range. Therefore, these two heat exchanges cancel each other and it is expected that after 300π min the system follows the same pattern as for the sinusoidal heat load. This behavior was indeed observed for all state variables. As example, we show the time evolution of the liquid phase volume in Figure 9.
An important aspect of this result is that the approach used here can handle discontinuities without losing consistency.

Example 5
Examples 1-4 are of differential index 1. In this example, we present a problem of differential index 2. In this case, the same system of Examples 2-4 is considered, but the pressure inside the tank is kept at 0.466 MPa, there is no heat load, and the evolution of vapor stream flowrate G is determined by solving the set of equations presented in Section 2. It should be pointed out that the vapor stream flowrate G is responsible for the higher differential index of the problem (index 2), due to the fact that G is a non-explicit algebraic variable. Figures 10 and 11 show the time evolution of G and T, respectively. In the first minutes, there is a greater accumulation of material inside the tank because the vapor Evolution of the liquid phase in the tank (Example 4). We consider two cases: a sinusoidal heat load Q · = sin (0.01t) (full line) and a piecewise heat load (dashed line), as defined in Table 4.
160π < t < 240π Q · = sin (0.01t) 320π < t < 400π Q · = sin (0.01t) stream flow is set initially equal to zero. This accumulation gives rise to an initial increase of temperature. Because of this, the vapor stream flowrate increases very rapidly to keep the pressure constant. This permits the temperature to decrease again until steady state is reached.

Example 6
In this last example, we present another problem of differential index 2. The same system of Examples 2-5 is considered, but the propane composition in vapor phase is maintained constant, equal to 0.193, and the inlet pressure is 0.65 MPa.
There is no heat load, and the evolution of the vapor stream flowrate G is determined by solving the set of equations presented in Section 2. Variable G is still responsible for the higher differential index of the problem (index 2). Example 6 differs from Example 5 due to the fact that, in Example 6, the additional condition has been applied to a differential variable (variable N 1 ), meanwhile in Example 5 the additional condition was applied to an algebraic variable (variable P 1 ). Figure 12 shows the time evolution of G. Figure 13 shows the time evolution of T and P.  Temperature inside the tank (Example 5).

Figure 13
Temperature and pressure inside the tank (Example 6).

CONCLUSIONS
The dynamics of separation drums was formulated as a system of differential algebraic equations (DAE), using the Peng-Robinson equation of state to calculate thermodynamic properties. We used PSIDE to integrate the DAE system, obtaining robust and efficient performance when solving the problems discussed here. This approach joins accumulated knowledge in numerical methods to rigorous thermodynamic treatment of the system. As demonstrated in Examples 5 and 6, we succeeded in the direct solution of problems with differential index higher than 1, which are usually a source of difficulty for typical software used to solve DAE systems.
Given the generality of the formulation, we believe that the procedure would work with other equations of state of similar complexity, such as many cubic EOS, and with other EOS that are non-cubic in molar volume, even though this has not been tested yet.