Phase stability analysis of liquid-liquid equilibrium with stochastic methods

THERMODYNAMICS AND SEPARATION PROCESSES

Phase stability analysis of liquid-liquid equilibrium with stochastic methods

G. Nagatani^I; J. Ferrari^I; L. Cardozo Filho^II; C. C. R. S. Rossi^II; R. Guirardello^III; J. Vladimir Oliveira^I; M. L. Corazza^I,*

^IDepartment of Food Engineering, URI, Campus de Erechim, Phone: +(55) (54) 3520-9000, Fax: +(55) (54) 3520-9090, Av. Sete de Setembro 1621, CEP: 99700-000, Erechim - RS, Brazil. E-mail address: [e-mail]

^IIDepartment of Chemical Engineering, Maringá State University, UEM, Av. Colombo 5790, CEP: 87020-900, Maringá - PR, Brazil

^IIIDPQ-FEQ, College of Chemical Engineering, State University of Campinas, UNICAMP CEP: 13081-970, Campinas - SP, Brazil

ABSTRACT

Minimization of Gibbs free energy using activity coefficient models and nonlinear equation solution techniques is commonly applied to phase stability problems. However, when conventional techniques, such as the Newton-Raphson method, are employed, serious convergence problems may arise. Due to the existence of multiple solutions, several problems can be found in modeling liquid-liquid equilibrium of multicomponent systems, which are highly dependent on the initial guess. In this work phase stability analysis of liquid-liquid equilibrium is investigated using the NRTL model. For this purpose, two distinct stochastic numerical algorithms are employed to minimize the tangent plane distance of Gibbs free energy: a subdivision algorithm that can find all roots of nonlinear equations for liquid-liquid stability analysis and the Simulated Annealing method. Results obtained in this work for the two stochastic algorithms are compared with those of the Interval Newton method from the literature. Several different binary and multicomponent systems from the literature were successfully investigated.

Keywords: Phase stability; Liquid-liquid equilibrium; Stochastic algorithms; Thermodynamic models.

INTRODUCTION

The phase split phenomenon in liquid mixtures is related to thermodynamic stability of the system and comprises important information for the projection, modeling and simulation of industrial processes that involve liquid phase separation. Phase stability analysis allows determination of the exact number of coexisting phases in stable equilibrium and also provides estimation of phase compositions, affording a suitable initialization for phase equilibrium calculations.

Phase stability analysis of multicomponent mixtures is usually investigated by the Gibbs tangent plane distance criterion (Baker et al., 1982; Michelsen, 1982). For a given temperature and pressure, the necessary and sufficient condition for a phase with composition z to be stable is that the Gibss free energy surface of the mixture is not intercepted by the tangent hyper plane associated with this surface in x. To fulfill this condition, the Gibbs tangent plane distance function, TPD(x), must be nonnegative for any acceptable x. Therefore, it is possible to examine whether a phase is stable by minimization of the TPD(x) function, subject to mass balance constraints. If the tangent plane distance function at the global minimum point has a value greater than or equal to zero, then the analyzed phase is thermodynamically stable because the TPD(x) function is also nonnegative for all x vectors contained in the permissible region. If the TPD(x) function is negative at its global minimum point, then the tangent plane lies above the Gibbs free energy surface, and hence the examined phase is unstable and will be subdivided into new phases.

For the computational evaluation of this test, two distinct approaches can be employed: (i) resolution of an algebraic nonlinear set of equations that represents the stationary points of the Gibbs tangent plane distance function and (ii) direct minimization of the TPD(x) function, taking into account the constraints imposed by the mass conservation principle. In both cases, the application of conventional mathematical models presents difficulties, such as dependence on an arbitrary initialization or the possibility that the solution converges to a trial root or to a local minimum. Since multiple roots may be found, especially when liquid phases are under inspection, the first strategy requires a mathematical method able to find all roots for the set of equations. The second approach, in contrast, requires a reliable global optimization method.

The majority of numerical methods applied to phase equilibrium calculations has local convergence characteristics resulting in only local solutions and are very sensitive to initial guesses. Some work is available in the current literature regarding application of deterministic techniques for global optimization, such as homotopy continuation, branch-and-bound and Interval Newton analysis (Sun and Seider, 1995; McDonald and Floudas, 1995; Souza et al., 2006; Lima et al., 2006). However, just a few studies make use of global stochastic techniques, like Simulated Annealing and Genetic Algorithm, for phase equilibrium and phase stability calculations (Zhu and Xu, 1999; Zhu et al., 2000; Rangaiah, 2001).

In this context, the aim of this work is to provide a systematic comparison between the two mentioned strategies for phase stability analysis: Simulated Annealing (Press et al., 1992), wich is a global optimization procedure, and a subdivision method (Smiley and Chun, 2001; Corazza et al., 2007) – henceforth denominated SubDivNL – employed for solution of the stationary points on TPD of liquid mixtures modeled by the NRTL activity coefficient model (Renon and Prausnitz, 1968). Comparison of the approach used in this work with Interval Newton analysis is also provided.

PHASE STABILITY ANALYSIS

The Gibbs free energy for liquid mixtures at low or moderate pressures can be modeled in terms of activity coefficients. Thus, the Gibbs tangent plane distance for a mixture with n number of components at specified temperature T and pressure P is given by (Baker et al., 1982):

where γ and are the activity coefficient of component i of tried and tested phases, respectively, and x = (x₁, ..., x_n) and z = (z₁, ..., z_n) represent the composition of tried and tested phases, respectively.

Since a phase is stable if TPD(x) > 0, the function represented in Eq. (1) must be minimized with respect to all possible compositions, i.e.:

with 0 < x_i< 1 and i = 1, ..., n.

As discussed previously, this optimization problem can be solved by direct minimization of the objective function (Eq. 1) or computation of the stationary points of the TPD(x) function, wich are determined by the first-order derivatives with respect to the (n – 1) independent mole fractions. The stationary points of TPD(x) can be obtained by solving the following set of algebraic nonlinear equations, subject to the mass balance restriction:

When the activity coefficients are estimated using the NRTL model, this set of equations may have multiple roots in addition to the trivial solution x = z.

Simulated Annealing

Simulated Annealing (SA) is a stochastic global optimization technique suitable for detecting a global minimum hidden among many local minima (Press et al., 1992). This method explores the analogy with the way that, when cooled slowly, a metal forms a crystalline structure with minimum energy. Metropolis et al. (1953) proposed that the thermal equilibrium of an annealing process could be simulated by the Boltzmann probability distribution.

The Simulated Annealing algorithm is based on the Monte Carlo method and makes use of an objective function φ(x) to replace energy E. The vector x is defined in the n-dimensional space of the optimization variables and represents a particular configuration of the system. A parameter T_A, analogous to the annealing temperature, is employed for evaluation of the Boltzmann factor. In accordance with Faber et al. (2005), a new vector, x^k+1, is randomly generated from the earlier state x^k. If φ(x^k+1) < φ(x^k), then the new vector is unconditionally accepted, i.e., x^k+1 = x^k = x^opt; otherwise, the acceptance probability is determined by the Metropolis criterion.

The progressive decrease in annealing temperature reduces the probability that a more energetic state be selected. As suggested by Faber et al. (2005), the annealing schedule can be set by for j = 1, ..., ζ, in which 0 < α < 1. However, just as excessively fast cooling may produce a structure with a higher energetic content, an inappropriate change in the annealing temperature can result in convergence to a local minimum. Furthermore, it is also essential to start from a conveniently high temperature in order to allow a more detailed exploration of the solution space of the problem.

Over the past years the SA algorithm has been applied to several types of phase equilibrium calculation problems involving phase stability tests and flash equilibrium calculations by minimizing the Gibbs free energy (Zhu and Xu, 1999; Zhu et al., 2000; Rangaiah, 2001; Teh and Rangaiah, 2003; Corazza et al., 2004; Souza et al., 2004; Bonilla-Petriciolet, 2006) and parameter estimation of thermodynamics models (Henderson et al. 2001; Moura et al., 2005; Carvalho Jr., 2006; Franceschi et al., 2006). Preliminary studies indicated that the best values for phase stability calculations are = 1000.0, α = 0.9 and z = 10000. In this work, the SA algorithm used was that presented in Press et al. (1992), where the probability acceptance criterion used is Pr(φ(x^K+1)) > xr, Pr is the Boltzmann probability distribution and xr is a random number [0,1].

Subdivision Algorithm

According to Smiley and Chun (2001) the subdivision algorithm consists in a given set of algebraic nonlinear equations F(x) and a primitive interval (hyper rectangle) R ∈ R^d, where d is the problem dimension, and the aim of a subdivision algorithm is to determine all values x* = {x ∈ R^d|F(x) = 0} through systematic subdivisions and trials in R. The first partition of R to congruent subrectangles is by division of all coordinates (variables) into two equal parts so as to obtain 2^d subrectangles R_ij from one rectangle in the R^d space (the index j identifies the subrectangles produced at a subdivision level i).

For each subdivision level, the subintervals obtained are tested for the existence of roots in order that only those containing one or more solutions are maintained. The selection criterion is based on calculation of the Euclidean norm of the set of equations:

where

then the subrectangle R_ij is retained or else, it is discarded (y is the vector of the randomly sampled points), and a_k and b_k are the lower and upper bonds of unknowns.

After a finite number of subdivisions, a conventional method (e.g. the Newton-Raphson method) is employed to determine the roots, taking as initial values the middle points of the last subintervals. Since the solutions are confined to the retained subintervals, the convergence of the adopted method becomes more efficient and reliable (Smiley and Chun, 2001).

A recent study, a subdivision algorithm for phase equilibrium calculations of systems at high-pressures modeled with an equation of state was proposed. The results obtained were quite satisfactory (Corazza et al., 2007). It may be interesting to note that, to our knowledge, an application of a subdivision algorithm to liquid-liquid equilibrium calculations is not available in the literature. In this work, from preliminary tests, the appropriate maximum number of subdivisions for phase stability tests was determined to be six.

RESULTS AND DISCUSSION

To validate the algorithms used in this work (Simulated Annealing and SubDivNL subdivision method), some results presented in the literature for liquid-liquid stability calculations with the NRTL model were selected. The investigated systems include binary, ternary and multicomponent mixtures. All the calculations were performed in a PC Pentium 4 2.66GHz with 512MB RAM.

Case 1: citric acid (1) + butan-2-ol (2) system

In Table 1 the results of phase stability tests for the citric acid (1) + butan-2-ol (2) binary system at 25ºC and different overall compositions are presented. The binary parameters of the NRTL model were taken from Gecegormez and Demirel (2005), which were presented by Lintomen et al. (2001). One can see in Table 1 that, for this system, solutions obtained in this work are in perfect agreement with the roots found by Gecegormez and Demirel (2005), which used an interval analysis algorithm. It can also be observed that the SubDivNL algorithm had lower computational effort in terms of CPU time that the SA algorithm (our results) and was shown to be much less time-consuming than the Interval Newton (IN/GB) algorithm reported on Gecegormez and Demirel (2005).

Case 2: n-pentanol (1) + 2,2-dimethylbutane (2) system

Results of phase stability analysis obtained in this work for the n-pentanol (1) + 2,2-dimethylbutane (2) system at 25 ºC and several global compositions are presented in Table 2. Results from the SA and SubDivNL algorithms are compared with solutions reported by Gecegormez and Demirel (2005). The NRTL parameters, which refer to a VLE set of parameters previously reported in the literature, were taken from Demirel and Gecegormez (1991).

It can be observed in Table 2 that for the two overall compositions, z = (0.05, 0.95) and z = (0.20, 0.80), results obtained in this work are different from those presented by Gecegormez and Demirel (2005). The TPD curves for mixtures containing respectively 0.05 and 0.20 in mole fraction of n-pentanol are represented in Figures 1 and 2. As clearly shown in these figures, the nontrivial solutions reported by Gecegormez and Demirel (2005) are not true solutions for the phase stability test, so the system has only one phase in both cases.

Moreover, the solutions presented in the literature are not stationary points on the TPD surface. In fact, the Euclidean norms of the set of equations calculated for such roots are far from zero: for z = (0.05, 0.95) and x = (0.1294, 0.8706), the value found is 0.0229, and for z = (0.20, 0.80) and x = (6.5086 × 10^-2, 0.9349), the result is 0.0292.

For all overall compositions tested for this system, our calculations only indicated the existence of a one-phase, homogeneous system, which just corroborates phase equilibrium data reported by Sayegh and Ratcliff (1976).

Nevertheless, one should call attention to the fact that the use of a set of parameters coming from VLE to LLE calculations is not recommended (Poling et al., 2000). In this work, however, a parameter set from VLE (Demirel and Gecegormez, 1991) for LLE computations is applied just to enable numerical comparison with the results of Gecegormez and Demirel (2005).

Case 3: water (1) + butyl glycol (2) system

The NRTL parameters for this system, which are presented in the literature (Demirel and Gecegormez, 1991) and refer to a set of parameters from VLE, were taken from Gecegormez and Demirel (2005). Solutions of the phase stability test in relation to LLE for the water (1) + butyl glycol (2) system at 5ºC and several global compositions are presented in Table 3. As depicted in Figures 3 and 4 the nontrivial roots reported by Gecegormez and Demirel (2005) for the water overall molar fractions of 0.05 and 0.10, respectively, are not true roots. According to the results from the SA and SubDivNL algorithms presented in this work, this system has only one stable phase at these compositions.

Case 4: water (1) + citric acid (2)+ 2-butanol (3) system

The LLE binary parameters of the NRTL model for this system were used as presented by Gecegormez and Demirel (2005), originally reported by Lintomen et al. (2001). Table 4 contains the results of phase stability analysis for this ternary system at 25ºC and different overall compositions. Solutions found in this work are compared with the roots presented by Gecegormez and Demirel (2005).

As shown in Table 4, some divergences between our results and those presented in the literature can be verified. For the global composition of z = (0.10, 0.05, 0.85), the composition with the lowest TPD value found in this work was x = (8.4640 × 10^-2, 0.4743, 0.4411) with TPD = -0.1609, while this root is not reported in the literature. Also, for the global composition z = (0.05, 0.10, 0.85) the SubDivNL and SA algorithms found an unstable composition, not reported in the literature. It can also be verified from Table 4 that for the last third overall composition, a second root, found in the literature, was not found by the SA and SubDivNL algorithms.

Another discrepancy was verified for the two overall compositions: a) z = (0.25, 0.05, 0.70) and b) z = (0.30, 0.05, 0.65) that correspond to mass fractions of a) z = (0.0682, 0.1456, 0.7862) and b) z = (0.0855, 0.1520, 0.7625), respectively. While unstable roots from the phase stability test are presented in the literature, indicating the formation of two liquid phases, application of the SA and SubDivNL algorithms did not result in the occurrence of separate liquid phases (see Table 4). In order to elucidate these results, the ternary experimental diagram for this system is presented in Figure 5. The experimental data were obtained from the literature (Lintomen et al., 2001), and as one can see, the two feed points (a) and (b) lie in a homogeneous region, which is in complete agreement with the results found in this work.

Case 5: acetonitrile (1) + benzene (2) + n-heptane (3) system

Table 5 contains the results of phase stability analysis for this ternary system at 45ºC and different overall compositions. Solutions found in this work are compared with the roots presented by Gecegormez and Demirel (2005). As shown in this table, application of the subdivision algorithm indicated the existence of a third root for all compositions analyzed. This additional root however is not presented in the work of Gecegormez and Demirel (2005). For this problem, the Simulated Annealing algorithm had a better performance in terms of CPU time than the subdivision method.

Case 6: n-propanol (1) + n-butanol (2) + benzene (3) + water (4) system

Solutions of the phase stability test for this system are presented in Table 6. The binary interaction parameters for the NRTL model were taken from Tessier et al. (2000). Aiming at reducing the computational effort, the maximum number of subdivisions was set to four. It can be verified that both algorithms employed in this work were able to find the same solutions as those presented in the literature (Tessier et al., 2000). The CPU time values listed in this table confirm the advantageous performance of the Simulated Annealing algorithm for systems with larger number of components.

Case 7: n-propanol (1) + n-butanol (2) + benzene (3) + ethanol (4) + water (5) system

Table 7 depicts the results of the thermodynamic stability test for this multicomponent system. For this system the NRTL parameters were also taken from Tessier et al. (2000), whose results are compared in this work. The maximum number of subdivisions was again set to four. Although the Simulated Annealing algorithm had a better performance in terms of CPU time than the SubDivNL, solutions obtained with both mathematical methods employed in this work are in agreement with those reported by Tessier et al. (2000), who used the Interval Newton method.

CONCLUSIONS

In this work thermodynamic phase stability with two stochastic global algorithms was examined for binary, ternary and multicomponent liquid mixtures using the NRTL activity coefficient model. Results of application of the subdivision algorithm to liquid-liquid equilibrium calculations show that it was able to find all roots of phase stability tests for the systems investigated with relatively short CPU times for binary and ternary systems. Simulated Annealing was shown to be reliable and robust for systems with larger numbers of components. In a general sense, the two algorithms are relatively simple, provide reliable solutions for phase stability analysis, are easy to implement and require relatively low computational effort.

NOMENCLATURE

d

problem dimension

(-)

f(x)

nonlinear equation

(-)

F(x)

system of nonlinear equations

(-)

R

multidimensional rectangle

(-)

a_k

lower coordinate of the "child" rectangles at "k" dimension

(-)

b_k

upper coordinate of the "child" rectangles at "k" dimension

(-)

A_k

lower coordinate of the "parent" rectangles at "k" dimension

(-)

B_k

upper coordinate of the "parent" rectangles at "k" dimension

(-)

NRTL

nonrandom two liquid

(-)

i

subdivision level

(-)

y

vector of random point in R_ij and vapor phase mole fraction

(-)

x _ij

middle point of the rectangle "j" at the subdivision level "i"'

(-)

x

vector of liquid phase mole fraction

(-)

z

vector of global mole composition

(-)

T

temperature

(-)

T₀

initial annealing temperature

(-)

iRan

sample number in R_ij

(-)

iMax

maximum subdivision number

(-)

iCov

maximum coverage number

(-)

TPD

tangent plane distance

(-)

Greek Letters

V_k

coordinate of the middle point of "parent" rectangle at "k" dimension

(-)

γ_i

activity coefficient of component "i"

(-)

τ_ij

subdivision parameter of the rectangle "j" at the subdivision level "i"

(-)

α_m,n

matrix of zero and one elements

(-)

Subscripts

0

pure component

(-)

n

number of components in the mixture and n^th components of the system

(-)

ACKNOWLEDGEMENTS

The authors thank CENPES/PETROBRAS, FAPERGS and CNPq for the scholarships and financial support of this work.

(Received: February 19, 2007 ; Accepted: March 5, 2008)

• Baker, L. E., Pierce, A. C. and Luks, K. D., Soc. Pet. Eng. J., 22, 731 (1982).
• Bonilla-Petriciolet, A., Vazquez-Roman, R., Iglesias-Silva, G. A. and Hall, K. R., Performance of stochastic global optimization in the calculation of phase stability analyses for nonreactive and reactive mixtures, Ind. Eng. Chem. Res., 45, 4764 (2006).
• Carvalho Jr., R. N., Corazza, M. L., Cardozo-Filho, L. and Meireles, M. A. A., Phase equilibrium for (camphor+CO₂), (camphor+propane), and (camphor +CO₂+propane), J. Chem. Eng. Data, 51, 997 (2006).
• Corazza, M. L., Corazza, F. C., Cardozo Filho, L. and Dariva, C., A subdivision algorithm for phase equilibrium calculations at high pressures, Braz. J. Chem. Eng, In press (2007).
• Corazza, M. L., Filho, L. C., Oliveira, J. V. and Dariva, C., A robust strategy for SLV equilibrium calculations at high pressures, Fluid Phase Equilibria, 221, 113 (2004).
• Demirel, Y. and Gecegormez, H., Simultaneous representation of excess enthalpy and vapor-liquid equilibrium data by NRTL and UNIQUAC models, Fluid Phase Equilibria, 65, 111 (1991).
• Faber, R., Jockenhövel, T. and Tsatsaronis, G., Dynamic optimization with Simulated Annealing, Computers & Chem. Eng., 29, 273 (2005).
• Franceschi, E., Kunita, M. H., Rubira, A. F., Muniz, E. C., Corazza, M. L., Oliveira, J. V. and Dariva, C., Phase behavior of binary and ternary systems involving carbon dioxide, propane, and glycidyl methacrylate at high pressure, J. Chem. Eng. Data, 51, 686 (2006).
• Gecegormez, H. and Demirel, Y., Phase stability analysis using Interval Newton method with NRTL model, Fluid Phase Equilibria, 237, 48 (2005).
• Henderson, N., Oliveira Jr., J. R., Amaral Souto, H. P. and Marques, R.P., Modeling and analysis of the isothermal flash problem and its calculation with the Simulated Annealing algorithm, Ind. Eng. Chem. Res., 40, 6028 (2001).
• Lima, E. R. A., Tavares, F. W., Biscaia, E. C., Jr., Utilização de Algoritmos Heurísticos de Otimização para Cálculo de Equilíbrio de Fases. In: XVI Congresso Brasileiro de Engenharia Química COBEQ, Santos (2006).
• Lintomen, L., Pinto, R. T. P., Batista, E., Meirelles, A. J. A. and Maciel, M. R. W., Liquid-liquid equilibrium of the water+citric acid+short chain alcohol+tricaprylin system at 298.15 K, J. Chem. Eng. Data, 46, 546 (2001).
• McDonald, C. and Floudas, C. A., Global optimization for the phase stability problem, AIChE J., 41, 1798 (1995).
• Metropolis, N., Rosenbluth, A. W. and Rosenbluth, M. N., Equation of state calculations by fast computing machines, J. Chem. Phys., 21, 1087 (1953).
• Michelsen, M., The isothermal flash problem. Part I. Stability, Fluid Phase Equilibria, 9, 1 (1982).
• Moura, L. S., Corazza, M. L., Cardozo-Filho, L. and Meireles, M. A. A., Phase equilibrium measurements for the system fennel (Foeniculum vulgare) extract + CO₂, J. Chem. Eng. Data, 50, 1657 (2005).
• Poling, B.E., Prausnitz, J.M. and O'Connell, J.P., The Properties of Gases and Liquids, 5^th ed., McGraw Hill, New York (2000).
• Press, W. H., Flannery, B. P., Teukolsky, S. A. and Veterlling, W. T., Numerical Recipes in FORTRAN - The Art of Scientific Computing, 2^nd ed., Cambridge University Press, New York (1992).
• Rangaiah, G. P., Evaluation of genetic algorithms and Simulated Annealing for phase equilibrium and stability problems, Fluid Phase Equilibria, 187-188, 83 (2001).
• Renon, H. and Praunitz, J. M., Local compositions in thermodynamic excess functions for liquid mixtures, AIChE J., 14, 135 (1968).
• Sayegh, S. G. and Ratcliff, G. A., Excess Gibbs energies of binary systems of isopentanol and n-pentanol with hexane isomers at 25 şC: Measurement and prediction by analytical group solution model, J. Chem. Eng. Data, 21, 1, 71 (1976).
• Smiley, M. W. and Chun, C. J., An algorithm for finding all solutions of a nonlinear system, Computational and Applied Math., 137, 293 (2001).
• Souza, A. T., Cardozo-Filho, L. and Wolff, F., Application of interval analysis for Gibbs and Helmholtz free energy global minimization in phase stability analysis, Braz. J. Chem. Eng., 23, 117 (2006).
• Souza, A. T., Corazza, M. L., Cardozo-Filho, L., Guirardello, R. and Meireles, M. A. A., Phase equilibrium measurements for the system clove (Eugenia caryophyllus) oil + CO₂, J. Chem. Eng. Data, 49, 352 (2004).
• Sun, A. C. and Seider, W. D., Homotopy-continuation method for stability analysis in the global minimization of the Gibbs free energy, Fluid Phase Equilibria, 103, 203 (1995).
• Teh, Y. S. and Rangaiah, G. P., Tabu search for global optimization of continuous functions with application to phase equilibrium calculations, Computers & Chem. Eng., 27, 1665 (2003).
• Tessier, S. R., Brennecke, J. F. and Stadtherr, M. A., Reliable phase stability analysis for excess Gibbs energy models, Chem. Eng. Science, 55, 1785 (2000).
• Zhu, Y. and Xu, Z., A reliable prediction of the global phase stability for liquid-liquid equilibrium through the Simulated Annealing algorithm: Application to NRTL and UNIQUAC equations, Fluid Phase Equilibria., 154, 55 (1999).
• Zhu, Y., Xu, Z. and Wen, H., Global stability analysis and phase equilibrium calculations at high pressures using the enhanced Simulated Annealing algorithm, Chem. Eng. Science, 55, 3451 (2000).

*

To whom correspondence should be addressed

Publication Dates

• Publication in this collection
  02 Sept 2008
• Date of issue
  Sept 2008

History

• Received
  19 Feb 2007
• Accepted
  05 Mar 2008