Local Thermodynamic Models Networks for Dynamic Process

Ind. Eng. Chem. Res. 2009, 48, 8529–8541

8529

Local Thermodynamic Models Networks for Dynamic Process Simulation Pedro Rafael Bolognese Fernandes* and Jorge Ota´vio Trierweiler Group of Integration, Modeling, Simulation and Optimization of Processes (GIMSCOP), Department of Chemical Engineering, Federal UniVersity of Rio Grande do Sul (UFRGS), Rua Luiz Englert s/n, CEP 90040-040, Porto Alegre, RS, Brazil

This paper presents the local thermodynamic models networks (LTMN), a new technique aimed at reducing the time spent in simulating processes that are based on phase equilibrium. The LTMN shown in this work approximate K-values as generated by a conventional property method, with a sensible reduction of the computer burden. The LTMN combine the local thermodynamic models technique and the concepts of multiple modeling to provide a valid approximation on the entire domain of the independent thermodynamic state variables. The paper discusses topics from the local models employed in the LTMN and the division of the space to the interpolation algorithm. Additionally, a new polynomial-time algorithm for the calculation of the bubble point temperature is developed on the basis of the form of the thermodynamic local models. The methods presented are tested in the dynamic simulation of distinct separation processes, for two distinct multicomponent mixtures, with time savings up to 90%. 1. Introduction The use of computer simulation as a process analysis tool is now well established in chemical engineering, ranging from design, static and dynamic simulation, to control analysis and optimization, among others, due to the valuable information it provides with reasonable realism at a low cost. The processes based on phase equilibrium, due to their economic significance, are regarded as an important field for the application of computer simulation techniques. A typical commercial simulation program for such processes usually contains a module, usually called as the “property package” or “physical property system”, destined to calculate the value of thermophysical properties (TP) that appear in the rigorous, phenomenological models. The description of phase equilibrium behavior by means of increasingly more accurate and sophisticated equations of state/activity coefficient models makes this property package complex, expensive, and computationally demanding. Moreover, the complex functional form presented by these thermodynamic models with respect to the state variables such as composition, temperature and pressure, often in implicit form, and the high number of gradient evaluations that most simulation algorithms need, frequently make the property package responsible for the largest contribution to the overall simulation time. In this way, Chimowitz et al.1 and Chimowitz et al.2 affirm that about 80% of the time spent in the dynamic simulation of phase equilibrium processes results from the calculation of thermophysical properties. Even with the availability of more powerful hardware, such kind of problems is still preventing the use of computer simulation in its full potential for industrial applications. Therefore, the reduction of the computational burden associated with the calculation of the TP can have a significant impact on the simulation time as a whole, and it can allow for applications that are otherwise impractical, as, for example, screening design, real-time optimization and online simulation (including model based control, real-time dynamic optimization, property estimation, etc.). Among the various strategies that have been proposed in the past for the reduction of the time spent in the calculation of * To whom correspondence must be addressed. E-mail: pedro@ enq.ufrgs.br.

thermophysical properties,3 the local thermodynamic models technique (LTM) is a straightforward and efficient alternative, and it is currently in use by several commercial process simulators (Aspen Engineering Suite, Hysys Process). The LTM appeared in the works of Hutchinson and Shewchuk,4 Leesley and Heyen,5 and Barret and Walsh.6 In addition, a closely related concept was successfully employed in a class of methods for the solution of the equilibrium stage problem, the inside-out algorithms,7 which employs local models for accelerating the convergence during the initial steps of the calculation. Basically, the LTM is a simple approximation function that is fitted by means of a few parameters to data generated with the conventional property model. Among the properties of interest, the Volatilization equilibrium ratio or (K-value), is essential for phase equilibrium based unit operations: Ki )

sat yi φˆ Li γiPsat i φi ) V ≈ xi φˆ i φˆ Vi P

(1)

where yi and xi are respectively the mole fractions of the component i in the vapor and liquid phases, γi is the liquid phase activity coefficient, Psat i the vapor pressure of component i pure, ˆV ˆL P is the system pressure and φsat i , φi and φi are the gas phase fugacity coefficients at the (pure) saturation condition and in the vapor and liquid phases of the multicomponent mixture, respectively. The LTM with best results were derived from phase equilibrium relationships;6,1,5 the simulation time reduction ranged from 30 to 60% in the literature, what may be interesting for some applications. In particular, one typical K-value local model, appropriate for the description of nonideal behavior due to the introduction of composition dependence through simplified activity coefficient terms (pseudobinary mixture hypothesis1,3) has the form ln(KiP) ) θi,1 +

θi,2 + θi,3(1 - xi)2 T

(2)

where T is the absolute temperature and θi,j, i ) 1, ..., Nc (number of components), j ) 1, 2, 3 are the adjustable model parameters. For fixed values of T and xi, and with Ki calculated by means

10.1021/ie9001469 CCC: $40.75  2009 American Chemical Society Published on Web 08/07/2009

8530

Ind. Eng. Chem. Res., Vol. 48, No. 18, 2009

of a conventional property method, it is possible to determine the parameters θi,j such that the local model (eq 2) is able to provide a good approximation to Ki, as long as the independent variables do not deviate too much from their previous values. In general, nevertheless, the periodic correction, or “update”, of these parameters is critical for the accuracy of the approximation since the model is inherently “local”, in the sense that it represents the phase equilibrium adequately only close to the region in which it was adjusted. Several strategies3,8,9 have been proposed in the literature to cope with the question of parameter correction during the simulation. However, some authors10,11 found problems with the use of the LTM with modern, multiple step DAE (differential-algebraic equation) system solvers because of the parameter discontinuity after the update. The approach discussed in this paper does not suffer from this problem, since the parameter discontinuity is eliminated, as discussed in the next section. The records from the published literature on the topic of approximate functions for speeding up dynamic simulations almost vanish by the 2000s, what may reveal a decreasing interest on analytical methods in face of increasing hardware capacities or, perhaps, that the subject is at a mature stage of development. Nonetheless, the development of more complex thermodynamic models and the online application of process simulation still justify the interest on the topic. 2. Local Thermodynamic Models Networks Although the LTM have been used in commercial process simulators with satisfactory results as an alternative to the calculations of TP with conventional models, this methodology has some problems: • Because of its local nature, parameters of the LTM have to be stored for each component in each equilibrium stage, thus increasing the total amount of information to be handled and complicating the implementation of the algorithm. • The LTM technique is not independent of the simulation program, since the solver has to “know”, among all the local models, which one has to be called for a given equilibrium stage in a given moment. • It is very difficult to determine a priori when the parameters should be updated, since the actual approximation error is not known, except when the conventional property routine is called. As a result, it is difficult to find the best compromise between precision and time reduction.3 In the same way, it is not easy to know what amount of past information has to be discarded for calculating the next set of parameter values. • The update procedure will in general cause discontinuities in the values predicted by the local models. This fact can be problematic when using the LTM in dynamic simulation.10,11 The local thermodynamic models networks (LTMN) technique was devised to overcome these inherent problems of the original LTM approach.12-14 This new methodology combines the LTM with the concept of multiple modeling based on the operating regime approach.15,16 Instead of dynamically tracking the system behavior through frequent parameter update, the LTMN is composed by local models representing the property of interest in specific subregions of the thermodynamic space, consequently forming a network. The models corresponding to each subregion are combined, or “blended” to provide a global property model, that is, an approximation valid in the entire thermodynamic space. The resulting network needs not to be

Figure 1. Local thermodynamic model network approximating the K-value in a hypothetical mixture at a fixed pressure.

updated after it was constructed and constitutes a “stand-alone” model for the property and the given multicomponent mixture. Figure 1 illustrates the case where a LTMN was constructed to approximate the K-value of a given component in a ternary mixture. In such case, three local models of the form (2) with distinct sets of parameters constitute the network for the entire thermodynamic space, here defined in terms of liquid composition. An interesting feature of the LTMN is that is represents the property for each component in the whole process flowsheet, unlike the LTM technique, in which one local model has to be stored for each different compartment, such as a single equilibrium stage of a distillation column. Furthermore, the LTMN, due to the interpolation procedure, generates continuous and smooth predictions, avoiding the problems related in the literature with DAE integrators. Moreover, the LTMN allows for the estimation of the approximation error, since its construction is based on the calculation of equilibrium data with the conventional property method covering the thermodynamic region of interest. In the same sense, the LTMN can be constructed with a given accuracy, a higher precision implying in a larger number local models and consequently in a more complex network. Finally, a single LTMN can be generated by means of conventional thermodynamic models of different forms or even with distinct sets of models parameters, for different regions of the thermodynamic space. Basically, the steps needed for the construction of the network are as follows: • Acquisition of equilibrium data, which is calculated with the conventional thermodynamic models. They can be regarded as the “training set” or the “experimental data” on which a model must be fitted. • The space division procedure, which includes the determination of LTM parameters for each subregion. • The combination of the local models in order to generate the global model. After the last step the network is ready for use in the process simulation. The final form of the LTMN is a routine which accepts as inputs state variables such as x, y, P, or T and returns the value of the property for the given input variables. In the sequel, only LTMN for the K-values will be considered, but the method is still applicable with other properties of interest, given that an appropriated local model is used to describe the corresponding property. 3. Local Models for the LTMN The LTM of the form given by eq 2 are considered in this work, since local models derived from rigorous phase equilib-


rium relationships ensure various desirable characteristics, as for example, a larger region of validity than it would be possible with the use of arbitrary functions. This larger validity region permits the construction of LTMN with a smaller number of subregions, what represents computational time savings during network generation and property estimation. Moreover, this physical background represents advantage over other possible approximation techniques, such as neural networks and direct multidimensional interpolation/look-up schemes. For instance, the physical background permits the LTMN to be extrapolated with more accuracy than the neural networks; it makes possible the model parameters to be related to physically meaningful variables; and it ensures less model complexity. In this way, another possibility is the consistent generation of derivatives of thermodynamic properties (e.g., ∂Ki/∂P, ∂Ki/∂T, etc.). Moreover, an interesting feature of model 2 is its linear structure in the parameters and in transformed independent variables, which is very amenable to the mathematical treatment of the problem. Additionally, the model represented by eq 2 does not require the calculation of other thermodynamic properties, like vapor pressure, for example, what is an advantageous characteristic. 3.1. New Thermodynamic Local Model. The thermodynamic local models considered so far in the literature were relatively simple functions, penalizing the number of parameters in order to require less data storage. The periodic update approach employed in these methods makes possible the satisfactory use of local models containing few parameters, since in general only a small region in the around the state trajectory direction will be spanned until the next update. Analyzing the functionality of eq 2, it must be noted that the model can be considered “symmetric”, in the sense that, for a fixed value of xi, any distribution of the remaining components of the mixture will result in the same value for the local model (pseudobinary mixture hypothesis). This means that the component i is assumed to interact in a similar way with all components of the mixture, what is, of course, not valid in general. Therefore, a new thermodynamic local model is proposed here, motivated by the fact that a better description of the composition effects can lead to a smaller number of subregions and consequently simplify the network. In this sense, the following local model for a multicomponent mixture is proposed: ln(KiP) ) θi,1 +

θi,2 + T

Nc

∑

θi,kxj2, k ) 3, 4, ..., Nc + 1

j)1,j*i

(3) where Nc is the number of components. This model has Nc 2 more coefficients than the local model in eq 2, therefore for binary mixtures the models are identical. The capability of the local models given by eqs 2 and 3 for modeling the phase equilibrium of some mixtures is further discussed in section 7. Another difference with other previous approaches is that no reference component6 was considered in this work, since it tends to make the algorithm more complex and leads to problems when the system composition is substantially changed.9 3.2. Adjusting the Local Models. The model parameters θi,j can be determined from a given set of Np equilibrium points calculated with the conventional routines trough an ordinary least-squares procedure if new independent variables are conveniently defined (for example, 1/Ti and (1 - xi,1)2 for model 2), and considering the transformed dependent variables ln(KiP) or ln(Ki). Optionally, one could use the best uniform approximation criterion, what can be accomplished through a linear

8531

programming problem. It should be nevertheless noted that the approximation (according to any norm) in terms of ln(KiP) does not necessarily give the same result as the approximation in terms of Ki. It was assumed here, however, that the loss incurred in using this representation is small in comparison with the augmented complexity of solving a nonlinear optimization problem. 4. Partitioning the Thermodynamic Space The fundamental idea behind the LTMN is the division of the thermodynamic space such that the underlying function can be adequately approximated in each subregion by means of a simple local model, and thus, the combination of the submodels produces a good global representation. The division procedure should therefore yield a good representation of the phase equilibrium surface with a small number of subregions (models). The space division is essentially a procedure that seeks for the reduction of the approximation error through the successive inclusion of new local models in the network. Thus, it is very important that enough information about the system be available during the process in the form of equilibrium data, to correctly drive the division; on the other hand, the collection of too much data can be excessively time-consuming. For simplicity, it was considered in this work that all data, obtained through bubble point calculations with conventional thermodynamic models, was available before the division procedure took place, although other strategies are also possible. 4.1. Selection of Independent Variables. The operating regime approach entails a previous selection of the variables that describe the operating space (that is, the independent variables). In the case under analysis, not all 2Nc + 2 intensive state variables P, T, and phase compositions (x and y, from which only Nc - 1 is truly independent) can be freely changed, but according to the Gibbs’ phase rule, only Nc, among them, is independent for a vapor-liquid equilibrium state. On the basis of the thermodynamic consistency and the functionality of local models 2 and 3, P and x were chosen as independent dimensions defining the division space, and T and y were considered dependent variables. The choice of P was merely a matter of convenience, since in most process simulations it is more natural to specify the pressure than the temperature. 4.2. Algorithms for the Space Division. Among the partitioning algorithms, the division of the space in terms of geometric subregions, mainly hyperquadrilaterals, is probably the simplest choice and numerous techniques have been proposed in literature, such as LOLIMOT (local linear model trees17). This kind of algorithm is supposed to generate a local models network with a prescribed accuracy trough the successive division of the space into hyper-rectangles and the collocation of new local models. In each step of the algorithm, LOLIMOT divides at the middle the current subregion (that with the largest approximation error) orthogonally with respect to an independent variable, and two new local models (one for each of the new resulting subregions) are fitted to the data contained in these sets. This procedure is done successively for all the independent variables, and the division that mostly reduced the global error is retained, thus increasing the overall number of subregions by one. A new subregion is then chosen, and the whole process is repeated until a prescribed approximation tolerance is reached. In this work, another algorithm was proposed: LOLIOPT (Local Linear Model Optimized Trees). Instead of dividing each subregion at the middle, LOLIOPT determines an optimal coordinate for this division. The objective of this procedure is to ensure a smaller number of subregions than it would be found

8532


by LOLIMOT. The following scheme outlines the two algorithms (the LOLIMOT algorithm was adapted here): (1) Select a subregion for the division (the i-th region). This region is initialized as the whole operating space, including all experimental data. (2) For each of the Nd dimensions of the i-th hyperquadrilateral, do the following steps, that is, for j ) 1, ..., Nd: (a) determine a division coordinate δi,j for the i-th hyperquadrilateral in the j-th dimension (b) divide the i-th hyperquadrilateral along the j-th dimension at the coordinate δi,j (c) adjust a LTM for each of the two new hyperquadrilaterals (d) compute the resulting global approximation error of the network (3) Retain the division that reduced most the global error. (4) Replace the current subregion with the two new regions. (5) Determine which subregion has the largest error; this will be the current division region (i-th hyperquadrilateral). (6) Return to step 1 while the error criterion is not satisfied. For the rigid partition approach adopted here, step 2d) above is equivalent to the choice of the partition that resulted in the smallest error for the current subregion, since the local division procedure has no effect on the rest of the network. This also avoids that local parameters are optimized with information concerning the global error, what is a more complex problem subject to local minima and so on. Moreover, for the application under study here the network was supposed to minimize the largest deviation between the data and the local model within the subregion. The difference between LOLIMOT and LOLIOPT resides solely in the step 2a, that is, the determination of the coordinate δ at which the direction is to be divided. In the case of LOLIMOT, for an operating space defined in terms of variables z ∈ RNd, δi,j ) (zi,jsup - zi,jinf)/2, where zi,jsup, zi,jinf are respectively the upper and lower limits of the i-th hyperquadrilateral in the j-th dimension. Alternatively, LOLIOPT determines the best division coordinate δ by means of an iterative procedure (line search). Recently, other optimality criteria have been presented, for example, on the basis of particle swarm optimization,18 which can yield more efficient networks. 4.3. The h-Transformation. As stated in section 4.1, the thermodynamic independent variables chosen to define the operating space were the composition x and the pressure P. Nevertheless, the composition domain is not suited for the division procedure in terms of hyperquadrilaterals because of its geometry (that is, not defined by a Cartesian product but by the manifold ΣNcxi ) 1 on RNc). An algebraic transformation, based on the h-transform or hyperplane transformation,19 was therefore utilized to obtain an orthogonal space analogous to the composition domain and suited for the division algorithm. The h-transform is usually employed in the study of nonlinear wave propagation phenomena in separation processes. Given a mixture with composition [x1, x2, ..., xNc], the corresponding h-coordinates are the Nc - 1 roots of the following equation (here slightly rearranged from the original form into a polynomial): Nc

Nc

∑ ∏ (h - R xi

i)1

m,k)

)0

(4)

m,m*i

where the constant Ri,k was considered to be the relative volatility Ki/Kk of the component i with respect to a fixed component k in the original paper. The inverse h-transformation can be found

by a similar expression. The process of projecting a point of the composition space on the h-space thus implies in finding the roots of a polynomial of degree Nc - 1. It must be noticed that the h-transform is used here just to generate a transformed composition domain that is more convenient to the numerical algorithm and not for linearizing the system behavior. In this sense, the system relative volatities Ki/Kk need not to be constant, as in the original paper, and the Ri,k are in fact rather arbitrary, and can be for example considered to be the average volatilities of the entire data set, or the volatilities of an equimolar mixture. 5. Generating the Network In a rigid partition approach, the continuity and smoothness of the network through the boundaries of each subregion cannot be guaranteed. For this approximation scheme, the result after the space division is a piecewise model of the form

{

Φ(z;θI), z ∈ RI Φ(z;θII), z ∈ RII f(z) ≈ l Φ(z;θNr), z ∈ RNr

(5)

where RI, RII, ..., design the Nr disjoint subregions determined by the space division algorithm, and θI, θII, ..., are the parameters that fit the local model to f(z) within each region. If one admits that the parameters, as given in the representation above, are in fact a function g, RNd f RNp, to the set of network inputs in the set of local models parameters, then a suitable interpolation function can be found in order to represent the parameter variation across the boundaries of the subregions f(z) ≈ Φ(z;g(z)) ≡ g(z)Tz

(6)

Furthermore, if g is continuous and smooth, it is easy to verify that the network will be continuous and smooth, if for example the local models are linear (as assumed in eq 6). Application to the LTMN. Following this methodology, a global model for the K-value can be provided if an appropriated function of the independent variables, representing the parameter values, is available ln(KiP) ) θ(h)Tφ

(7)

with the variables of the thermodynamic local model included in the vector φ ≡ [1, 1/T, x12, xj2, ..., xNc2] (j * i) for model given by eq 3, and the model parameters θ ≡ [θi,1, θi,2, ..., θi,Nc+1] represented as a function of h, the vector of h-coordinates. The values of the parameters for each subregion, resulting from the division procedure, are assigned to the geometric center of each hyperrectangle, and a multidimensional interpolation procedure is applied to find a functional relationship between the parameter and the coordinates of the space. Moreover, the parameters θi,j, j ) 1, 2, ..., Nc, for a specific component i have the same network configuration, which is advantageous for the interpolation procedure. It must be noted that after the space division and the determination of the parameters of the interpolation functions, the network may be written in the form of an independent numerical routine that takes as inputs T, P and x and returns the corresponding value of K, thus making it suitable for computer simulation. 5.1. Modified Shepard Interpolation Scheme. Shepard interpolation scheme20-22 was considered in this work due to its favorable characteristics of simplicity and capability of


handling scattered data as generated by the partitioning algorithm. The basic idea of the Shepard method is that the value of an unknown function f at a given point z ) [z1, z2, ..., zNd]T can be approximated as the weighted sum of the known function values at the nodes ci, where ci ) [ci,1, ci,2, ..., ci,Nd]T, i ) 1, 2, ..., Np Np

f(z) ≈

∑ w (z)f(c ) i

i

(8)

i)1

where wi(z) are the weight functions. Particularly, the “localtype” version (with Franke-Little weights) is very suited to the type of application intended here, since it considers only data within a prescribed radius to calculate the weighted interpolation. To avoid the problems of defining the local radii, from which the interpolation quality is very dependent, and to take advantage of all information about the subregions obtained through the division procedure, a modification of the Shepard method is proposed here. In this case, the weights are given by wi(z) )

Φi(z)Ψi(z) Np

∏φ

j i,k(zk)

(12)

k)1

where k indexes the dimension, and

{

1 - 6ξji,k5 + 15ξji,k4 - 10ξji,k3 0 e ξji,k e 1

φji,k(ξji,k)

) 0

ξji,k > 1

1

ξji,k < 0 (13)

The transformed position variable ξ is given by ξji,k(zk) )

zk - ck,i Fji,k - ck,i

(14)

where Fjj,k is the maximum coordinate of influence of the region j into region i in the k-th dimension, which is determined by 2-β (ck,j - Rji,k) + Rji,k 2

(15)

(9) j represents the coordinates of the boundary between and Ri,k regions i and region j in the direction k, given by

i

i)1

where Ψi are the nodal functions, that is, functions that quantify the influence of the nodal points on the weighted interpolation. In this work, an exponential function of the following form was considered: Ψi(z) ) exp(-γ·di(z)2)

φji(z) )

Fji,k )

∑ Φ (z)Ψ (z) i

8533

Nd

(10)

The scalar factor γ is a tuning parameter of the nodal function, and the function di is the -norm of the vector z - ci, defined in a generalized manner to represent the rectangular shape of the subregions. The Φi are the internodal functions that quantify the allowable influence of a given node on the others and are employed to ensure a smooth transition between the subregions. These functions have many aspects in common with the blending functions described by Mo¨llers et al.,23 although they were considered here in a more multidimensional setting. The modification presented in this work was called complete blending function because the influence of all subregions is considered to compute the weight associated with a given node. Therefore, it is necessary to consider blending functions associated with a given pair i-j of subregions; the cumulative effect of this blending on region i from the adjacent regions is given by

Rji,k ) ck,j + lj sign(ck,i - ck,j)

(16)

which is specified the network configuration. The smoothness of the function is controlled by means of the scalar β e 2. Figure 2 exhibits an example of interpolation in one dimension using this method. A set of evaluations of a given function f(x), x ⊂ R, at the centers corresponding to 4 subregions was interpolated with the scheme above. The example shows the main advantage of interpolation functions of this class, namely, that the resulting function is smooth and monotone, which is highly desirable when interpolating the parameters of the local models at the centers of the subregions. 6. Bubble Point Temperature Approximation The functional form of eqs 2 and 3 permits one to obtain expressions that approximate the bubble point temperature explicitly, what avoids the need of solving a set of nonlinear equations through a recursive algorithm to determine the phase equilibrium. If the thermodynamic local model 2 is placed in the “bubble point equation”

Nneig{i}

Φi(z) )

∏

φji(z)

(11)

j)1,j*i

To guarantee the local character of the interpolation scheme, this product is done over the number of neighboring regions of i, Nneig{i}, defined as the number of regions whose boundaries have at least one point in common with the boundary of region i. Since these blending functions must consider the relative position of each subregion, the resulting functions are in general multivariable. The pairwise blending function associating the influence of region j in the region i, φij(z), is the product of unidimensional functions

Figure 2. Modified Shepard’s interpolation scheme: interpolation of a set of points with two different β-γ parameter values.

8534

Ind. Eng. Chem. Res., Vol. 48, No. 18, 2009 Nc

∑ K (P, T, x , y )x i

i

i

i

)1

(17)

i)1

∑

1 Kixi ) P i)1

∑ exp(θ Nc

i,1

+

i)1

Na

∑ j!1 q (T - T )

j

j

0

(23)

j)1

where T0 is a convenient reference temperature (for example, (Ta + Tb)/2), and the coefficient qj is given by

the result is Nc

( )

( )

θl,2 θl,2 ≈ exp + exp T T0

)

θi,2 + θi,3(1 - xi)2 · xi ) 1 T (18)

( )∑

θl,2 · qj ) exp T0

(θl,2)i j!(j - 1)! (-1) i!(i - 1)!(j - i)! Ti+j i)1 0 j

[

j

]

(24) Local model in eq 2 was employed above for shortness; similar results will apply to model in eq 3. The solution of eq 18, given the pressure and the liquid phase composition, is the bubble point temperature of the system, approximated by means of a thermodynamic local model. If we assume that the term exp(θi,2/T) can be represented by an explicit function of the temperature, for example

exp

( )| θi,2 T

T)Tb T)Ta

≈ bi,0 + bi,1T

(19)

then the bubble point temperature will be directly given by

7. Example: Construction of a Network

Nc

P-

∑

bi,0 exp(θi,1 + θi,3(1 - xi)2)xi

i)1

T≈

(20)

Nc

∑b

i,1

exp(θi,1 + θi,3(1 - xi)2)xi

i)1

For better accuracy, the approximation above was made in a suitable interval [Ta, Tb], which can correspond to the minimum and maximum pure-liquid boiling points at the pressure P. In the case of systems forming azeotropes, information about the maximum/minimum equilibrium temperature can be used or, in the absence of that, a convenient range around [Ta, Tb]. The accuracy of this approximation or, equivalently, the linearity of exp(θi,2/T), will depend on θi,2; for this reason, higher approximation orders are generally necessary. The result is then a polynomial expression of the type Na

Nc

∑ ∑b Tk

k)0

i,k

A schematic representation of the bubble point temperature approximation algorithm is shown in Figure 3. It should be noted that, although the least-squares method gives the best least-squares fit for low approximation orders (1-6), this kind of approximation can be very unstable as the order increases, depending the length of the interval [Ta, Tb] and on the value of θi,2. Consequently, some kind of control of the approximation order must be kept to avoid such problems. In opposition, the use of Taylor series permits the selection of arbitrarily high approximation orders, but in this case, the higher computational load should also be considered.

exp(θi,1 + θi,3(1 - xi)2)xi ) P

(21)

i)1

where Na is the approximation order. One of the roots of this polynomial constitutes an approximate solution T* of the original bubble point problem; this root must be physically meaningful and satisfy, to a given accuracy, the summation of mole fractions Nc

∑ K (T*)x i

i

≈1

(22)

i)1

The approximation coefficients bi,k are not additional parameters of the LTMN technique because they are determined exclusively by the values of θi,j that result from the network (represented by functions of the h-coordinates). Several methods can be used to determine the bi,k values: for example, the expansion of exp(θi,2/T) in Taylor series, the polynomial leastsquares fit (both continuous and discrete), and so on. Since the function to be approximated is analytically known, it is possible to obtain the approximations in closed form in both cases. For example, the Taylor approximation of exp(θl,2/T) around T0 is

The methodology described in the previous sections and all necessary thermodynamic routines were implemented in Matlab, and a LTMN was constructed to approximate the K-values of the system water, acetonitrile, and acetone at pressures around 1 atm. This system is considerably nonideal and was employed earlier in the literature on local thermodynamic models8,9,11 as a benchmark example. For this purpose, 441 bubble point problems, corresponding to a uniform grid in the h-space, were calculated. The h-coordinates were defined in the range 1 e h1 e 1.284 and 1.284 e h2 e 2.686, which are the relative volatilities of the components of this mixture with respect to water in an equimolar mixture at 1 atm, as discussed in Section 4.3. The conventional thermodynamic models were constituted by the Soave-Redlich-Kwong equation of state for the vapor phase and by the Wilson model with temperature-dependent parameters for the liquid phase. Vapor pressure was calculated either with the Frost-Kalkwarf-Thodos (acetonitrile) or the Wagner (water and acetone) equations. The parameters used in the conventional thermodynamic routines can be found in Reid et al.24 For each of the components, LOLIOPT was employed to divide the h-space until an approximation tolerance of 3% in the K-values was reached. This deviation refers to the maximum distance between the data and the adjusted piecewise local models and does not take into account the effect of the interpolation procedure. Local models of form in eq 3 were utilized, and the best uniform approximation was employed to fit local models parameters to the data, resulting in 15 subregions for water, 5 for acetonitrile, and 4 for acetone. The least-squares criterion yield a 50% larger number of centers. The network configuration resulting of this procedure is given in Figure 4. When compared to LOLIOPT, LOLIMOT did not result in a remarkably worse number of centers (16, 7, and 7). Moreover, the use of model in eq 2 instead of model in eq 3 gave as result 6 (water), 21 (acetonitrile), and 20 (acetone) centers. The comparison of the prediction capability of models in eqs 2 and 3 in this and in some other systems (Table 1) shows that model 3 generally outperforms model in eq 2. Nevertheless, the network configuration for the water was evidently worse for model in eq 3. This fact suggests that the K-value of water in this system depends much more on the water composition than


8535

Figure 3. Polynomial-time bubble point temperature approximation algorithm.

Figure 4. Result of the space division in h-coordinates for the ternary system at 1 atm: water (left), acetonitrile (center), and acetone (right). Table 1. Maximum Relative Deviation in K (%) when Adapting Models 2 and 3 to 441 Equilibrium Points at 1 atm mixture

component

model 2

model 3

1

n-pentane n-hexane n-heptane acetone benzene ethanol methyl ethyl ketone ethanol water

2.17 3.30 3.37 12.03 100.23 65.60 127.13 57.74 63.71

0.69 0.41 2.75 28.99 31.54 25.38 49.95 31.00 61.86

2 3

on the relative amount of other components. This example illustrates nicely the fact that a model like eq 3, even with more degrees of freedom, can produce worse predictions if it does not correspond to the physics of the system. Analyzing Figure 4, it can be seen that due to the lack of prediction capability of model in eq 3, the LOLIOPT algorithm divided the space into very narrow regions with respect to the h2 axis. This led to a situation in which very few data points were contained in these subregions (i.e., as many points as model parameters), and as a result, the local models could not be fitted properly. In view of this, local eq 3 was replaced by eq 2 for water in this system. This also exemplifies one characteristic of the LTMN, namely, the possibility of using local models of different structures in the same network. Besides the time necessary to evaluate the equilibrium data, LOLIOPT needed about 1 min per component to generate the network in this ternary example (that is, a bidimensional search in h-space), while LOLIMOT took less than 10 s per component. For a quaternary example, the time spent by LOLIOPT raised to 10 min, whereas LOLIMOT need about 30 s per component. These local models were then employed to generate a LTMN for the system water, acetonitrile and acetone by means of the method proposed in section 5. The modified Shepards interpolation method was used, with interpolation parameters β ) 0 and γ ) 0.1, which are supposed to give a smooth resulting network.

7.1. Bubble Point Calculations. This LTMN was then compared with the conventional thermodynamic model in the computation of 961 bubble points. For this data set, the mean relative deviations in the K-values were 2.2% (water), 1.2% (acetonitrile), and 0.72% (acetone), while the maximum relative deviations were 6.7%, 7.4%, and 12.2%, respectively. The result of the bubble-point calculations are shown in Figure 5. The maximum deviations in the K-values are larger than the specified tolerance (3%) because in the division procedure the temperature is given (as an independent variable), and in the case of bubble point calculations the temperature is determined on the basis of the K-value model, thus propagating the error of the network. The mean and the maximum deviation in the bubble point temperature were 0.35 and 1.8 K, respectively. For this kind of calculation, the LTMN model allowed a reduction of about 90% in the computation time when compared to the conventional model (from 928 to 75 s for the 961 points). If the Frost-Kalkwarf-Thodos equation, which is implicit, is replaced by the Antoine equation for the acetonitrile in the conventional model, this saving is reduced to 35% since the vapor pressure calculation is very time-consuming with the first method. Furthermore, the ideal model Ki ) Pisat/P (Raoult’s law, with explicit pressure vapor equations) gives essentially the same computational gain as the LTMN models, but in this case, the maximum and mean deviation in temperature were 15.8 and 6.6 K, respectively. 8. Application to Process Simulation The local thermodynamic models networks generated for two ternary systems, namely, water-acetonitrile-acetone and benzene-ethanol-acetone, were employed in the simulation of some simple separation problems (dynamic flash vessels and distillation column with equimolar flows). Simulations were carried out in a PC with AMD Athlon 64 1.8 GHz processor and 450 Mb of RAM.

8536


Figure 5. Top: K-values of water with conventional thermodynamic models (left) and LTMN (right). Center: K-values of acetone with conventional thermodynamic models (left) and LTMN (right). Bottom: Bubble point temperature with conventional thermodynamic models (left) and LTMN (right), for the water-acetonitrile-acetone system at 1 atm.

8.1. DAE Flash Problem. System Water, Acetonitrile, and Acetone. The LTMN constructed in the previous section for the system water, acetonitrile and acetone was tested in the dynamic simulation of a flash vessel modeled as a DAE system (the model is given in the Appendix). Vessel pressure was controlled at a constant reference value of 1 atm by manipulating the vapor outflow, and the feed stream was a mixture with composition zf at its bubble point at 1.5 atm. For this test, 30 different runs were considered, consisting in a step change in feed composition from zf(a) to zf(b), with both zf(a) and zf(b) randomly determined in Matlab in the range 0.05 e zf,i e 0.95, i ) 1, 2. The dynamic response of the flash vessel was then simulated with the conventional routines and the LTMN. Integration was done by means of a variable order solver based on the numerical differentiation formulas (Matlab v. 5.3/6.1 function ode15s25) with absolute and relative tolerances of 10-3. Consistent initial conditions were generated by the solution of an equivalent stationary adiabatic flash problem with P ) 1 atm and feed composition zf(a), using the conventional thermodynamic models (for all simulations). It should be noted that this initial state introduces an extra difficulty for the initialization of the DAE problem with the LTMN, since small inconsistencies appear. In despite of that, the average computational gain of the LTMN over the conventional routines was 93.6%. Average simulation times were 8.25s (conventional method) and 0.48s

(LTMN). The results of two typical runs (zf(a) ) [0.1 0.1 0.8] and zf(b) ) [0.1 0.8 0.1], zf(a) ) [0.1 0.1 0.8] zf(b) ) [0.8 0.1 0.1]) and are shown in Figure 6. System Benzene, Ethanol, and Acetone. To check the results, the benzene, ethanol and acetone system was taken as example. The thermodynamic space was divided with the LOLIOPT algorithm by means of 441 equilibrium points at 1 atm obtained with the Soave-Redlich-Kwong EOS, the UNIQUAC activity coefficient model and explicit vapor pressure equations, with a tolerance of 2%, giving a number of centers of [5 2 4]. A LTMN was constructed with the interpolation scheme discussed in section 5. The results of the network for the prediction of the equilibrium temperature are given in the Figure 7. For comparison, the result obtained with Raoult’s law is also presented. For the bubble point problem, the average computational gain of the LTMN with the multidimensional blending functions was 55%. The ideal model needed 44% less time then the conventional model in this case. The use of this LTMN with the DAE flash simulation tests showed an average gain of 25%; the average simulation times were 0.5s (rigorous model) and 0.33s (LTMN). 8.2. Simplified Flash Problem. The dynamic simulation of a simplified flash problem was used to check the behavior of the LTMN with a constant step ODE integrator (explicit Runge-Kutta 2,3 method implemented in the Matlab vs 5.3/ 6.1 function ode23, absolute and relative tolerances of 10-3).


8537

Figure 6. Dynamic response of the flash vessel modeled as a DAE system with respect to two step changes in feed composition: liquid mole fraction (top), vapor mole fraction (center), and temperature (bottom). Reponses with the conventional model are shown by the continuous lines, and responses with the LTMN are shown by the dotted lines.

The model is given in the Appendix. Two step changes in feed composition of the benzene, ethanol and acetone system were simulated: from [0.1 0.8 0.1] to [0.8 0.1 0.1] and from [0.1 0.8. 0.1] to [0.1 0.1. 0.8] (Figure 8). The time reduction for the LTMN was in the order of 40-45%. The average simulation time was 2.1s for the conventional method and 0.85s for the LTMN. 8.3. Bubble Point Approximation Algorithm. Equilibrium Calculations. The bubble point approximation method of section

6 was tested with this system. The approximation interval was chosen to be [329.2, 353.3] K, which are the pure boiling points of acetone and benzene at 1 atm, calculated with the Wagner equation. On the basis of the LTMN for this system, the analysis of the maximum absolute values of the parameter θi,2 among all subregions (approximately -2.8 × 104 for benzene, -3.5 × 103 for ethanol, and -9 × 103 for acetone), suggested a Taylor expansion of order 7. The bubble point algorithm with Taylor approximation required 77% less time than the bubble

Figure 7. Bubble point temperature for the mixture benzene, ethanol and acetone at 1 atm: conventional routine (SRK-UNIQUAC, left), LTMN (center), and ideal (right).

8538


Figure 8. Responses of liquid-phase mole fraction in face of two feed changes, according to the conventional thermodynamic model (continuous lines) and to the LTMN (dotted lines).

Figure 9. Bubble point temperature of the system benzene, ethanol, and acetone at 1 atm with the approximation algorithm (left); liquid phase mole fraction response of the simple flash problem according to the conventional thermodynamic model (continuous lines) and the LTMN with the bubble point temperature algorithm (dotted lines).

point problem with the LTMN, with small temperature deviations. When compared to the problem with conventional thermodynamic methods, the time-savings reached 90%. Figure 9 shows the result of the Taylor approximation to the bubble point temperature. Simplified Flash Problem. The approximation algorithm with Taylor series was also tested in the dynamic simulation of the simplified flash problem; the computational time gain in calculating the response of a step change in feed composition (from [0.1, 0.8, 0.1] to [0.1, 0.1, 0.8]) was 75%, corresponding to average simulation times of 2.1s (conventional) and 0.5s (LTMN with bubble point algorithm). This result can be seen in Figure 9. Dynamic Simulation of a Distillation Column. The LTMN for the system benzene-ethanol-acetone with the Taylor approximation algorithm was also tested in the dynamic simulation of a distillation column with equimolar flows, that is, with vapor phase composition and equilibrium temperature calculated by means of the bubble point problem. The model is

described in the Appendix. The column is constituted by 20 trays, and a 10 mol/s saturated liquid mixture of benzene, ethanol, and acetone at 1.5 atm is fed in the stage 4 from the top. A linear pressure profile was considered, with condenser pressure of 1 atm. For this test, a change of feed composition, from zf (a) ) [0.3 0.2 0.5] to zf (b) ) [0.5 0.2 0.3], was simulated with the conventional routines and the LTMN. Integration was done by means of the Matlab vs 5.3/6.1 function ode23tb (implicit Runge-Kutta formula) with absolute and relative tolerances of 10-3. The computational gain of the LTMN over the conventional routines was around 90%, with simulation times of 250s (conventional method) and 25s (LTMN). The result of the feed composition step change can be seen in Figure 10. 9. Discussion and Conclusions This paper presented a new technique destined to the reduction of the computational burden associated with the

Figure 10. Dynamic profile responses of the distillation column model: stage temperature (left), benzene vapor phase mole fraction (center), and acetone liquid phase mole fraction (right), according to the conventional thermodynamic model (continuous lines) and to the LTMN (dotted lines).


simulation of liquid-vapor equilibrium separation processes. In despite of the increasing hardware speed, dynamic simulation of these separation processes is still a challenge, in the sense that many important tasks (parameter estimation, simulation of control strategies, online optimization, and so on) are still prohibitive for most applications. Moreover, modern equations of state based on descriptions at molecular or quasi-molecular level, as for example the SAFT and COSMO-RS types, are still very computer intensive for engineering applications. The examples showed that the LTMN approach can lead to significant time savings which range from 40% to 95%, depending on the nature of the problem (static/dynamic), the formulation of the model (DAE or ODE), the nature of mixture under study, the form of the network, among others. In particular, it was found that the LTMN suffers to less extent from the problems found in the literature for the simulation of DAE systems with local models techniques, which led to an increase of the solution time. Nevertheless, an increase in the number of integration steps with the LTMN was found in all simulations, and also a small increase in the number of failed steps. The basis of the LTMN technique is that a global model of the property can be constructed trough the combination of several, piecewise defined, local models. This is obviously an empirical assumption, but it can be expected to work in practical situations such as in the reconstruction of continuous, smooth phase equilibrium surfaces. Naturally, phenomena such as the formation of a second liquid phase, for example, cannot be predicted, unless the conventional model takes it into account and the network is precise enough. For simplicity, it was considered in this work that all data was obtained previously to the space division procedure, what is not the unique possibility; for example, “adaptive” division schemes can de devised, in which the division algorithm itself determines the region to be explored. This kind of approach can demand less data, since it will not be uniformly distributed, and regions with more complex phase equilibrium behavior will be reproduced in more detail. A problem with this approach is that no initial error estimate is available to drive the mechanism of placing the centers. Another possibility is the construction of the network in parallel to the simulation, as the thermodynamic space is being covered with data; after a certain period, the LTMN would be adequately constructed and could then be used for the calculation of thermophysical properties in substitution to the conventional model. One of the intrinsic problems of the type of approximation method pursued here is the question of dimensionality, often called “curse of dimensionality”. The amount of data needed to cover uniformly the operating space grows exponentially with the number of dimensions. Some schemes may be devised to alleviate this problem for systems with many components, as for example techniques for dimensionality reduction (SVM, PCA, PLS), use of coarser grids for less important components, grouping components, etc. The LTMN were constructed and simulated, for simplicity, at a fixed pressure. Some tests nevertheless showed that the quality of the approximation can decay significantly away from this reference pressure. As already stated in the literature,18 a solution could be the inclusion of an additional term of the form θp · (P/Pref) in the local model. Appendix. DAE Flash. The model of the DAE flash vessel is given by the Nc + 1 differential equations

8539

dNi ) F · zi - L · xi - V · yi dt

(A.1)

dU ) F · hf - L · h - V · H + Q dt

(A.2)

where Ni is the total number of moles of component i, F is the feed rate (10 mol/s), zi is the feed composition, L and V are the outlet liquid and vapor mole flow rates (mol/s), U is the total internal energy (J), hf, h, and H are the enthalpies of feed, liquid, and vapor stream, respectively, and Q is the heat added (considered zero in all simulations). The 2 Nc + 2 algebraic equations are given by yi - Kixi ) 0 Nc

∑

(A.3)

Nc

yi -

i)1

∑x

i

)0

(A.4)

i)1

Nc

Ni - ni -

∑ (N

i

- ni) · yi ) 0

(A.5)

i)1

Nc

U-

∑ (N

Nc

i

- ni) · H -

i)1

∑n

i

· h - 0.0242PVvessel ) 0

i)1

(A.6) where ni is the number of moles of species i in the liquid phase and Vvessel is the volume of the flash vessel (6000 L). It was considered also that the outflows are controlled according to the equations V ) Vbias + 10 · (P - 1 atm) and L ) Lbias + 10-3 · (Vl - 3000 L), where Vbias and Lbias are the vapor and liquid flows at the stationary operating condition, P is the actual pressure of the vessel (calculated by an EOS as a function of the state variables), and Vl is the actual liquid phase volume (calculated with the COSTALD correlation20). K-values and enthalpies are calculated with a thermodynamic routine (conventional or LTMN). Local models for enthalpies are considered to be the enthalpies or an ideal solution. State variables were scaled to improve the numerical integration. Simplified Flash Problem. The model of the simplified flash vessel is given by the following Nc differential equations: dxi 1 ) (F · zi - L · xi - V · yi), i ) 1, ..., Nc - 1 dt n (A.7) dxi )dt

Nc-1

∑ j)1

dxj , i ) Nc dt

(A.8)

where n is the total holdup in the liquid phase (assumed to be constant), the pressure is assumed to be perfectly controlled at 1 atm (V is constant), L is given as in the DAE model, and the yi values are calculated for each integration step through a bubble point problem as a function of xi and P ) 1 atm, with the corresponding thermodynamic routine (conventional/LTMN/ ideal). Simplified Distillation Column Model. The model of the distillation column considers equimolar flows, Ns equilibrium stages 100% efficient, negligible vapor holdup, no heat losses, and constant condenser/reboiler holdups. Stage temperature and vapor compositions are calculated by means of a bubble point problem. Equations:

8540


liquid phase molar holdups, j-th stage, except feed tray (Ns · Nc - 3 equations) dni,j ) Vj+1 · yj+1 + Lj · xj - Vj · yj - Lj · xj dt

(A.9)

liquid phase molar holdups, feed tray (Nc equations) dni,j ) F · zi + Vj+1 · yj+1 + Lj · xj - Vj · yj - Lj · xj dt (A.10) condenser D ) V2 - R

(A.11)

B ) LNs-1 - W

(A.12)

Vj ) W

(A.13)

reboiler

internal flows

nj - nj,base β

Lj ) Lj,base +

(A.14)

constitutive relations xi,j ) ni,j /

∑n

Literature Cited i,j

(A.15)

i

nj )

∑n

i,j

y ) mole fraction of the component i in the vapor phase of a mixture z ) feed composition z ) point in R z ) point in the RNd space zi,jsup,zi,jinf ) upper/lower coordinates limiting the i-th hyperquadrilateral in the j-th dimension Φ ) internodal function of the interpolation algorithm Φ(z;θ) ) (linear) multivariable function of z with parameters θ Ψ ) nodal function of the interpolation algorithm γ ) parameter of the nodal function Ri,k ) relative volatility of component i with respect to component k β ) lower bound of the normalized interval δi,j ) division coordinate in the j-th dimension for the i-th hyperquadrilateral (LOLIMOT/LOLIOPT algorithms) φisat ) fugacity coefficient of pure gas at the saturation condition φˆ iV, φˆ iL ) fugacity coefficients of component i in the vapor and liquid phases of the multicomponent mixture γi ) liquid phase activity coefficient of component i φ ) vector of variables of the thermodynamic local model Fi,kj ) parameter of influence of region i into region j in the k-th dimension θi,j ) j-th parameter of the local model for component i ζi,kj ) transformed distance variable of region i into region j in the k-th dimension

(A.16)

i

where ni is the number of moles of species i in the liquid phase, L (V) is the mole flow of liquid (vapor), and D (B) is the distillate rate (bottoms rate). Nomenclature bi,k ) k-th parameter for the component i in the approximation of the exponential function (bubble-point algorithm) ci ) center/node of the i-th subregion (interpolation point) d ) the norm function of the distance vector z - c f(z) ) scalar-valued multidimensional function f(z) ) scalar-valued unidimensional function g ) vector-valued multidimensional function h, hi, h ) h-transform variable, i-th coordinate in h-space, vector of h-coordinates Ki ) volatilization equilibrium ratio Nc ) number of components in a mixture Nd ) number of dimensions in the partitioning algorithm Nneig{i} ) number of neighboring regions of i Np ) number of equilibrium points for constructing the LTMN Nr ) number of disjoint subregions P ) pressure Pisat ) vapor pressure of component i pure Rk ) k-th subregion determined by the space division algorithm Ri,kj ) coordinate of the boundary between regions i and j in the k-th direction R ) the real set T ) temperature Ta, Tb ) limits of the interval for the bubble-point temperature approximation algorithm w ) weight function of the interpolation algorithm x ) mole fraction of the component i in the liquid phase of a mixture

(1) Chimowitz, E. H.; Anderson, T. F.; Macchietto, S.; Stutzman, L. F. Local models for representing phase equilibria in multicomponent, nonideal vapor-liquid and liquid-liquid systems. 1. Thermodynamic approximation functions. Ind. Eng. Chem. Proc. Des. DeV. 1983, 22, 217. (2) Chimowitz, E. H.; Anderson, T. F.; Macchietto, S.; Stutzman, L. F. Local models for representing phase equilibria in multicomponent, nonideal vapor-liquid and liquid-liquid systems. 2. Application to process design. Ind. Eng. Chem. Proc. Des. DeV. 1984, 23, 609. (3) Hillestad, M.; Sorlie, C.; Anderson, T. F.; Olsen, L.; Hertzberg, T. On estimating the error of local thermodynamic modelssA general approach. Comput. Chem. Eng. 1989, 13, 789. (4) Hutchison, H. P.; Shewchuk, C. F. A computational method for multiple distillation towers. Trans. Inst. Chem. Eng. 1974, 52, 325. (5) Leesley, M. E.; Heyen, G. The dynamic approximation method of handling vapor-liquid equilibrium data in computer calculations for chemical processes. Comput. Chem. Eng. 1977, 1, 109. (6) Barret, A.; Walsh, J. J. Improved chemical process simulation using local thermodynamic approximations. Comput. Chem. Eng. 1979, 3, 397. (7) Boston, J. F.; Sullivan, S. L. An improved algorithm for solving the mass balance equations in multistage separation processes. Can. J. Chem. Eng. 1972, 50, 663. (8) Macchietto, S.; Chimowitz, E. H.; Anderson, T. F.; Stutzman, L. F. Local models for representing phase equilibria in multicomponent, nonideal vapor-liquid and liquid-liquid systems. 3. Parameter estimation and update. Ind. Eng. Chem. Proc. Des. DeV. 1986, 25, 674. (9) Perregaard, J. Model simplification and reduction for simulation and optimization of chemical processes. Comput. Chem. Eng. 1993, 17, 465. (10) Ledent, T.; Heyen, G. Dynamic approximation of thermodynamic properties by means of local models. Comput. Chem. Eng. 1994, 18, Suppl. S87S91. (11) Støren, S.; Hertzberg, T. Local thermodynamic models applied in dynamic process simulation. Trans. Inst. Chem. Eng. 1994, 72, 395. (12) Bolognese Fernandes, P. R. Redes de Modelos Termodinaˆmicos Locais. MS Dissertation, UFRGS, Porto Alegre, Brazil, 2001. (13) Bolognese Fernandes, P. R. ; Trierweiler, J. O.; Secchi, A. R. Local Thermodynamic Models Networks: A Novel Approach for Process Simulation. Preprints of DYCOPS-6; Cheju Island, South Korea; 2001; pp 579-584. (14) Bolognese Fernandes, P. R. ; Trierweiler, J. O.; Secchi, A. R.; Wada, K. Local Thermodynamic Models Networks for Process Simulation. Preprints of EMPROMER 2001, Santa Fe´, Argentina; 2001; Vol. 1, pp 595600. (15) Johansen, T. A.; Murray-Smith, R. Multiple Model Approaches to Modelling and Control; Taylor & Francis: London, 1997.

Ind. Eng. Chem. Res., Vol. 48, No. 18, 2009 (16) Trierweiler, J. O.; Neumann, U. Local Models Networks: a simple solution for complex problems. Preprints of of COBEQ’98; Porto Alegre, Brazil; 1998. (17) Nelles, O. LOLIMOT- Lokale, Lineare modelle zur identifikation nicht-linearer, dynamischer systeme. Automatisierungstechnik 1997, 163. (18) Mehran, R.; Fatehi, A.; Lucas, C. Araabi, B. N. Particle Swarm Extension to LOLIMOT. Preprints of the Sixth International Conference on Intelligent Systems Design and Applications; Jinan, China; 2006; Vol. 2, pp 969-974. (19) Ting, J.; Helfferich, F. G.; Hwang, Y.; Graham, G. K. Experimental study of wave propagation dynamics of multicomponent distillation columns. Ind. Eng. Chem. Res. 1999, 38, 3588. (20) Basso, K.; Zingano, P. R.; Freitas, C. M. Interpolation of scattered data: Investigating alternatives for the modified Shepard method. Preprints of of SIBGRAPI’99; Campinas, Brazil; 1999. (21) Engeln-Mu¨llges, G.; Uhlig, F. Numerical Algorithms with C; Springer-Verlag: Berlin, 1996.

8541

(22) Shen, C. Y.; Reed, H. L.; Foley, T. Shepard’s interpolation scheme for solution-adaptive methods. J. Comput. Phys 1993, 106, 52. (23) Mo¨llers, T.; Laak, O. On the global approximation and interpolation of locally described real valued functions. Appl. Math.Comput. 1998, 93, 1. (24) Reid, R. C.; Prausnitz, J. M. ; Poling, B. E. The Properties of Gases and Liquids, 3rd ed.; McGraw-Hill: New York, 1987. (25) Shampine, L. F.; Reichelt, M. W. The Matlab ODE Suite. SIAM J. Sci. Comp. 1997, 18 (1), 1.

ReceiVed for reView January 27, 2009 ReVised manuscript receiVed June 11, 2009 Accepted July 7, 2009 IE9001469

Local Thermodynamic Models Networks for Dynamic Process

Recommend Documents