QUADRATURE ALGORITHMS FOR PHASE EQUILIBRIUM OF CONTINUOUS MIXTURES

Abstract

Several methods for computing the Gauss-Christoffel quadrature used for the adaptive characterization of continuous mixtures were compared as to their efficiency and robustness. Two mixtures with molar fraction distribution given by truncated gamma distributions were used. We analyzed the Product-Difference, the Golub-Welsch, the Long Quotient-Modified Difference and the Chebyshev algorithms using regular and generalized moments, when applicable. The robustness and computational efficiency of changes in the distribution variable and in the orthogonal polynomial family used to calculate the generalized moments were analyzed. The methods using generalized moments proved to be more robust than those that use regular moments. Although they are computationally more expensive, this cost increase is just around 20% for the Chebyshev algorithm. The resulting adaptive characterization was employed to solve the adiabatic vapor-liquid flash of these mixtures. The results showed that eight pseudocomponents were able to well represent the properties of the equilibrium streams, showing the high accuracy of this method.

Keywords:
Continuous thermodynamics; Complex mixtures; QMoM; Gauss-Christoffel quadrature; Phase equilibrium

INTRODUCTION

Several chemical processes deal with mixtures whose components cannot be fully identified, such as petroleum and polymeric solutions. This difficulty is related to the amount of components in the mixtures and to the proximity of their chemical and physical properties. According to Briesen and Marquardt (2003), a petroleum mixture can contain over 10⁶ components, making it impossible to characterize all these components.

Therefore, approximate techniques have been developed to characterize these mixtures and to solve the thermodynamic equations related to the phase equilibria. The most common form of dealing with this problem is to use a molar fraction distribution function associated with at least one characterization variable and then perform the calculations. For hydrocarbon mixtures, especially for homologous series, the molar mass is often enough to characterize the mixture properties (Lage, 2007). Other variables such as normal boiling temperature or the API degree of the mixture are also used for this purpose (Whitson and Brule, 2000).

The molar fraction distribution functions of continuous mixtures are usually described by well-known distribution functions from the literature, such as gamma, beta or exponential distributions (Huang and Radosz, 1991; Whitson and Brule, 2000). Then, the parameters of these distributions are estimated in order to adjust the probability density distribution function to the mixture properties. Experimental analysis techniques, such as chromatography and true boiling point curves, are used to achieve these estimations (Whitson and Brule, 2000).

Once the molar fraction distribution function is characterized, the thermodynamic equations can be used to solve engineering problems. However, due to the non-linearities of the thermodynamic models, the cases for which there is an analytic solution are rare. One example is the solution of an isothermic flash for a continuous mixture considering ideal solutions and ideal gas behavior (Lage, 2007). For general cases, the distribution function must be discretized into pseudocomponents before the thermodynamic equations can be numerically solved. Several methods for its discretization have been studied in the literature. The simplest methods are those that employ an equally spaced or a function-based discretization (Huang and Radosz, 1991), which are very simple to implement. However, the properties of the continuous mixture are accurately approximated by those of the discrete mixture only if a large number of pseudocomponents is used with these methods. Huang and Radosz (1991) showed that gamma and truncated Gaussian distributions need at least 40 pseudocomponents when using a uniform discretization, or 20 pseudocomponents when a logarithmic function was used for this purpose. Nevertheless, due to the peculiarities of each mixture, the choice of a function that best describes its changes is not straightforward, making this method not very efficient.

Gauss quadratures are found in the literature as more efficient alternatives to the above describe methods (Cotterman and Prausnitsz, 1985; Haynes and Matthews, 1991; Liu and Wong, 1997). As most distributions are semi-infinite, the Gauss-Laguerre quadrature is usually the natural choice. This methodology allows the use of less discretization points than in the previous method and proved to be very efficient when characterizing distributions. However, this method has fixed quadrature points, not being suitable to characterize changes in the distribution function, either in time or space (Briesen and Marquardt, 2003).

In order to characterize mixtures whose distribution may vary in time or space, Lage (2007) developed a characterization method for continuous mixtures based on the Gauss-Christoffel quadrature using the moments of the molar fraction distribution function. This method is analogous to the QMoM (Quadrature Method of Moments), developed by Mcgraw (1997) for the solution of aerosol dispersions. Applying the QMoM for the continuous thermodynamics equations, it allows the streams to be characterized by their moments and, whenever these moments change, the method calculates the new Gauss-Christoffel quadrature.

Lage (2007) showed that this method presents good results for the solution of the flash equations and showed how to deal with the mixing of two streams based on their moments. Petitfrere et al. (2014) tested the method described by Lage (2007) for the solution of the flash equations of real mixtures, proving its adequacy and accuracy. Extending the application of this method to distillation columns, Rodrigues et al. (2012) developed a method for the sequential simulation of a distillation column for separation of continuous mixtures. Whenever a stream suffers any kind of change during the simulation, the method was able to generate a new quadrature rule that adaptively characterizes this stream using modified pseudocomponents.

The core of the QMoM calculation is a routine to compute the Gauss-Christoffel quadrature. In the development of the QMoM for the solution of the populational balance equations, Mcgraw (1997) used the Product Difference Algorithm (PDA) developed by Gordon (1968). However, in the literature, several methods for the generation of this quadrature can be found. Lage (2007) and Rodrigues et al. (2012) used an implementation of the PDA where critical computations were performed in higher precision by separating the floating point numbers into their mantissa and exponent (Lage, 2007). This method showed to be more robust than the PDA developed by Gordon (1968). Petitfrere et al. (2014) used the Chebyshev algorithm (Chebyshev, 1858), implemented by Gautschi (1994) in the library ORTHPOL, concluding that this algorithm is more robust than the original PDA (Gordon, 1968).

John and Thein (2012) numerically analyzed several algorithms for computing the Gauss-Christoffel quadrature. They compared the PDA from Gordon (1968), the Golub-Welsch Algorithm (GWA), developed by Golub and Welsch (1969), and the Long Quotient-Modified Difference Algorithm (LQMDA), developed by Sack and Donovan (1972). The LQMDA is one of the two existing variants of the Chebyshev algorithm (Chebyshev, 1858). The other variant is the Wheeler algorithm (Wheeler, 1974). John and Thein (2012) concluded that the LQMDA (or any of its variants) is more robust and more efficient than the other two methods.

However, John and Thein (2012) used the LQMDA only with regular moments, even though this method (and its variants) supports the usage of generalized moments. This type of moment is more flexible to characterize continuous mixtures than the regular moments, which raises the expectation that the methods using them can be more efficient and robust. Therefore, in this work we extend the analyses carried out by John and Thein (2012) including the generalized moments and the effect of the orthogonal polynomial family used to generate the quadrature. The method was applied to the adiabatic flash equations in order to analyze the efficiency of the quadrature points to represent the molar fraction distribution function for the streams in equilibrium.

METHODOLOGY

The methods analyzed in this work were the PDA (Gordon, 1968), the GWA (Golub and Welsch, 1969), the LQMDA (Sack and Donovan, 1972) and the Chebyshev algorithm (Chebyshev, 1858), which can be found in Upadhyay (2012). All of them were applied using regular moments and the last two were also tested using generalized moments.

In order to analyze these methods that generate the Gauss-Christoffel quadrature, a computational program was implemented that, given a molar fraction distribution function, calculates the moments of this distribution, generates its quadrature rule and, when required, solves the adiabatic flash equations for the composition of the pseudocomponents of the equilibrium streams. The procedures and equations used in this program are detailed in this section. The development of these methods can be found in the original references and also in John and Thein (2012). Therefore, their details are not given here.

Moments calculation

Consider a molar fraction distribution, f(x), in terms of a characterization variable, x, defined in the interval [a, b]. This distribution can be characterized by its regular moments, defined as:

$${\mu }_k=\underset{a}{\overset{b}{{\displaystyle \int }}}{x}^{k}f\left(x\right)dx,k=\mathrm{0,1,2,...}$$ (1)

where k is the order of the moment. These moments carry information of the distribution that is used by the algorithms to generate the desired Gauss-Christoffel quadrature.

There are also the generalized moments of the distribution that are obtained by substituting the monomial x ^k in Equation (1) by the k-order polynomial, p _k (x), of an orthogonal polynomial family, resulting in:

$${\mu }_k^{\left(P\right)}=\underset{a}{\overset{b}{{\displaystyle \int }}}{p}_k\left(x\right)f\left(x\right)dx,k=\mathrm{0,1,2,...}$$ (2)

where µ_k ^(P) is the generalized moment of order k associated with this orthogonal polynomial family, which can be obtained by the following recursion formula:

$$\begin{array}{l}{p}_{-1}\left(x\right)=0\hfill \\ p_0\left(x\right)=1\hfill \\ x{p}_k\left(x\right)=a_k{p}_{k+1}\left(x\right)+b_k{p}_k\left(x\right)+c_k{p}_{k-1}\left(x\right),k=\mathrm{0,1,...}\hfill \end{array}$$ (3)

where the coefficients a _k , b _k and c _k define the given family. It is important to note that the generalized moments are reduced to the regular moments for a _k = 1 and b _k = c _k = 0.

Once calculated the n quadrature points, the first 2n regular or generalized moments from the distribution can be reconstructed respectively from:

$${\mu }_k={\displaystyle \sum }_{i=1}^n{x}_i^{k}f\left(x_i\right),k=\mathrm{0,...,2}n-1$$ (4)

$${\mu }_k^{\left(P\right)}={\displaystyle \sum }_{i=1}^n{p}_k\left(x_i\right)f\left(x_i\right),k=\mathrm{0,...,2}n-1$$ (5)

Computation of the quadrature

The studied algorithms generate the quadrature [x _j , ω_j ]ⁿ _j=1 , that optimally approximates the integral:

$$\underset{a}{\overset{b}{{\displaystyle \int }}}g\left(x\right)f\left(x\right)dx\approx {\displaystyle \sum }_{j=1}^{n}g\left(x_j\right){\omega }_j$$ (6)

where f(x) is the molar fraction distribution function and g(x) is a given property of the mixture component x with infinitesimal molar fraction f(x)dx. By interpreting Equation (6) as a discretization of the mixture, the abscissas, x _j , represent pseudocomponents and the weights, ω_j , are their molar fraction.

The inner product between two functions p(x) and q(x) with respect to the weight function f(x) in the interval [a, b] is defined as:

$$\langle p,q\rangle \equiv \underset{a}{\overset{b}{{\displaystyle \int }}}p\left(x\right)q\left(x\right)f\left(x\right)dx$$ (7)

As for any Gaussian quadrature rule, the abscissa of the quadrature given by Equation (6) are the roots of the n-order polynomial of the Christoffel orthogonal polynomial family using the inner product defined by Equation (7). However, since f(x) is the unknown mixture molar fraction distribution, it is impossible to have a priori knowledge of the recursion coefficients of the corresponding orthogonal polynomial family. Therefore, the quadrature rule must be numerically computed from the known moments of the distribution function, as detailed in the following.

These Christoffel orthogonal polynomials are generated by the following recursion:

$$\begin{array}{l}{P}_{-1}\left(x\right)=0\hfill \\ P_0\left(x\right)=1\hfill \\ P_{k+1}\left(x\right)=\left(x-{\theta }_k\right)P_k\left(x\right)-{\eta }_k^2{P}_{k-1}\left(x\right),k=\mathrm{0,1,...}\hfill \end{array}$$ (8)

where P _k (x) is the Christoffel polynomial of order k and the coefficients _k and θ_k can be obtained by:

$${\theta }_k=\frac{\langle x{P}_k,P_k\rangle }{\langle P_k,P_k\rangle },k\ge 0$$ (9)

$${\eta }_k^2=\left\{\begin{array}{ll}1\hfill & ,k=0\hfill \\ \frac{\langle P_k,P_k\rangle }{\langle P_{k-1},P_{k-1}\rangle }\hfill & ,k\ge 1\hfill \end{array}\right.$$ (10)

Combining the recursion equations from P ₀ to P _n−1 , the following system of linear equations can be assembled:

$$\left({\tilde{A}}_n-xI\right)P=b$$ (11)

where:

$${\tilde{A}}_n=\left(\begin{array}{llllll}{\theta }_0\hfill & 1\hfill & 0\hfill & \dots \hfill & \dots \hfill & 0\hfill \\ {\eta }_1^2\hfill & {\theta }_1\hfill & 1\hfill & \hfill & \hfill & ⋮\hfill \\ 0\hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & \hfill & ⋮\hfill \\ ⋮\hfill & \hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & 0\hfill \\ ⋮\hfill & \hfill & \hfill & {\eta }_{n-2}^2\hfill & {\theta }_{n-2}\hfill & 1\hfill \\ 0\hfill & \dots \hfill & \dots \hfill & 0\hfill & {\eta }_{n-1}^2\hfill & {\theta }_{n-1}\hfill \end{array}\right)$$ (12)

$$P=\left(\begin{array}{l}{P}_0\left(x\right)\hfill \\ P_1\left(x\right)\hfill \\ ⋮\hfill \\ P_{n-1}\left(x\right)\hfill \end{array}\right)$$ (13)

$$b=\left(\begin{array}{l}0\hfill \\ ⋮\hfill \\ 0\hfill \\ -P_n\left(x\right)\hfill \end{array}\right)$$ (14)

If x _i is one of the n roots of P _n , then the linear system reduces to an eigenvalue problem:

$$\left({\tilde{A}}_n-x_{i}I\right)P=0$$ (15)

Therefore, the process of finding the n roots of the orthogonal polynomial P _n (x) can be replaced by the process of finding the n eigenvalues of the matrix à _n . The numerical methods for finding the eigenvalues are better conditioned that those for finding the roots of polynomials, which represents a numerical advantage in exchanging these problems (John and Thein, 2012).

Furthermore, the computing of eigenvalues for symmetric matrices is more accurate than for non-symmetric matrices. Thus, the matrix à _n should be modified in order to create a symmetric matrix. Considering the diagonal matrix

$$D={[d_i]}_{i=0}^{n-1},$$

where

$$d_i={\left[{\displaystyle {\prod }_{j=0}^i{\eta }_j}\right]}^{-1}$$

for i = 0, …, n - 1, then:

$$D\left({\tilde{A}}_n-x_{i}I\right)D^{-1}DP=0$$ (16)

$$\left(\underset{A_n}{\underbrace{D{\tilde{A}}_n{D}^{-1}}}-x_i\underset{I}{\underbrace{DI{D}^{-1}}}\right)\underset{\widehat{P}}{\underbrace{DP}}=0$$ (17)

$$\left(A_n-x_{i}I\right)\widehat{P}=0$$ (18)

where

$$A_n=\left(\begin{array}{llllll}{\theta }_0\hfill & {\eta }_1\hfill & 0\hfill & \dots \hfill & \dots \hfill & 0\hfill \\ {\eta }_1\hfill & {\theta }_1\hfill & {\eta }_2\hfill & \hfill & \hfill & ⋮\hfill \\ 0\hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & \hfill & ⋮\hfill \\ ⋮\hfill & \hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & 0\hfill \\ ⋮\hfill & \hfill & \hfill & {\eta }_{n-2}\hfill & {\theta }_{n-2}\hfill & {\eta }_{n-1}\hfill \\ 0\hfill & \dots \hfill & \dots \hfill & 0\hfill & {\eta }_{n-1}\hfill & {\theta }_{n-1}\hfill \end{array}\right)$$ (19)

The eigenvalues of A _n , also called the terminal matrix, are the roots of the Christoffel orthogonal polynomial of order n, which are the abscissas of the Gauss-Christoffel quadrature.

After computing the eigenvector, q _i , associated with the eigenvalue x _i , the weights of the Gauss-Christoffel quadrature can be expressed by (John and Thein, 2012):

$${\omega }_i={\mu }_0^{\left(P\right)}{q}_{i\mathrm{,0}}^2{\left({\displaystyle \sum }_{k=0}^{n-1}{q}_{i,k}^2\right)}^{-1},i=\mathrm{1,...,}n$$ (20)

The sum inside Equation (20) is the norm of the eigenvector q _i . If the routine used to compute the eigenvectors is such that the eigenvectors are orthonormal, then this sum term can be removed from the equation.

The algorithms that compute the Gauss-Christoffel quadrature differ in how they build the terminal matrix. Once this matrix is obtained, the eigenvalue problem is solved in the same way for all methods. All methods use the moments of the distribution as the necessary information to build the matrix A _n . For a quadrature with n points, the first 2n moments must be known.

Considering the analyzed methods, it is important to highlight that only the LQMDA and the Chebyshev Algorithm can use both regular and generalized moments as input. Also, the PDA was tested in two different implementations: the original version by Gordon (1968) and the modified implementation given by Lage (2007).

Flash Equations

The equations related to the flash of a continuous mixture require the usage of continuous properties to solve the equations of mass balance, energy balance and phase equilibrium for the entire distributions. These equations are, respectively:

$$f^F\left(I\right)=\gamma f^V\left(I\right)+\left(1-\gamma \right)f^L\left(I\right),\forall I\in \left[I_{min},I_{max}\right]$$ (21)

$$H^F\left(T\right)=\gamma H^V\left(T\right)+\left(1-\gamma \right)H^L\left(T\right),\forall I\in \left[I_{min},I_{max}\right]$$ (22)

$$f^V\left(I\right)=K\left(I,T,P\right)f^L\left(I\right),\forall I\in \left[I_{min},I_{max}\right]$$ (23)

where I is the characterization variable, usually related to molar mass, and γ is the vaporized fraction. Assuming ideal solution and ideal gas, the enthalpy from the streams and the equilibrium constant can be calculated by Equations (24) and (25), respectively.

$$H^G\left(T\right)=\underset{I_{min}}{\overset{I_{max}}{{\displaystyle \int }}}{f}^G\left(I\right)H^{gi}\left(I,T\right)dI,G=F,L,V$$ (24)

$$K\left(I,T,P\right)=\frac{P^{sat}\left(I,T\right)}{P},\forall I\in \left[I_{min},I_{max}\right]$$ (25)

Once the Gauss-Christoffel quadrature is found for f ^F (I), the abscissas are the values of the characterization variable for the discretized pseudocomponents of the feed stream and the weights are their molar fractions. Then, the flash equations for discrete mixtures can be applied, keeping these pseudocomponents fixed during the flash calculation. Afterwards, the moment set of each resulting stream, liquid and vapor, can be used to generate a new pseudocomponent discretization by computing the corresponding Gauss-Christoffel quadrature.

The equations of an adiabatic flash of a discrete mixture consist of the mass and energy balances and the phase equilibrium equation, which are expressed, respectively, by Equations (26), (27) and (28).

$$x_i^F=\gamma x_i^V+\left(1-\gamma \right)x_i^L,i=\mathrm{1,...,}n$$ (26)

$$H^F\left(T\right)=\gamma H^V\left(T\right)+\left(1-\gamma \right)H^L\left(T\right),i=\mathrm{1,...,}n$$ (27)

$$x_i^V=K_i\left(T,P\right)x_i^L,i=\mathrm{1,...,}n$$ (28)

The equilibrium constants, K _i , and the enthalpy of the streams are calculated by:

$$K_i\left(T,P\right)=\frac{P_i^{sat}\left(T\right)}{P},i=\mathrm{1,...,}n$$ (29)

$$H^G\left(T\right)={\displaystyle \sum }_{i=1}^n{x}_i^G{\widehat{H}}_i^{G,id},G=F,L,V$$ (30)

where the ideal enthalpy for each component is given by:

$${\widehat{H}}_i^{V,id}=h_{f,i}^{gi}\left(T_0\right)+\underset{T_0}{\overset{T}{{\displaystyle \int }}}{C}_P^{gi}dT$$ (31)

$${\widehat{H}}_i^{L,id}=h_{f,i}^{gi}\left(T_0\right)+\underset{T_0}{\overset{T}{{\displaystyle \int }}}{C}_{P,i}^{gi}\left(T\right)dT-\Delta H_{vap,i}$$ (32)

where ΔH _vap,i was assumed to be independent of T and equal to its value at 298.15 K.

The empirical correlations used to estimate the thermodynamics properties were taken from the literature. These correlations relate the property from the pure component i to its molar mass, M _i . These correlations can be found in the Appendix or in the works of Huang and Radosz (1991) (saturation pressure) and Marano and Holder (1997) (enthalpy of formation, enthalpy of vaporization at 298.15 K and ideal gas heat capacity).

NUMERICAL PROCEDURE

Distribution Functions

For testing the different algorithms for Gauss-Christoffel quadrature computation, a molar fraction distribution is needed. Some petroleum mixtures have distributions that can be approximated by gamma distributions. Therefore, in this work, two truncated gamma distributions were used.

The gamma distribution is written as:

$$f\left(M\right)=\frac{1}{FN}\frac{{(M-M_0)}^{A-1}}{B^A\Gamma \left(A\right)}\text{exp}\left(-\frac{M-M_0}{B}\right)$$ (33)

where the parameters A, B, and M ₀ and the truncation interval for the two distribution functions are given in Table 1. These distributions are depicted in Figure 1. FN is a normalization factor due to the truncation, given by:

$$\text{FN}=\underset{\Delta M}{{\displaystyle \int }}\frac{{(M-M_0)}^{A-1}}{B^A\Gamma \left(A\right)}\text{exp}\left(-\frac{M-M_0}{B}\right)dM$$ (34)

Table 1
Parameters and truncation intervals for the two gamma distributions.

Distribution 1	Distribution 2
A = 2.1	A = 4.0
B = 26.7	B = 35.0
M0 = 100	M0 = 100
M = [100, 300]	M = [100, 450]

Figure 1
Distribution functions 1 and 2 used to analyze the quadrature algorithms.

Definition of the distribution variable

The distribution variable is defined to be a function of the molar mass. This variable change is important to increase the accuracy in the calculation of the moments and, therefore, the overall algorithm robustness. Lage (2007) commented that this change of variable is necessary to keep the moments with similar order of magnitude. This is done by defining I in order that it belongs to the [0, C] interval, where C is around 1.

Therefore, the distribution variable was defined by:

$$I=\frac{M-M_0}{M_f-M_0}C$$ (35)

and the corresponding molar fraction distribution was derived as:

$$f\left(I\right)=\overline{C}\frac{I^{A-1}}{{\overline{B}}^A\Gamma \left(A\right)}\text{exp}\left(-\frac{I}{\overline{B}}\right)$$ (36)

where

$$\overline{B}=\frac{BC}{M_f-M_0},\overline{C}=\frac{1}{\text{FN}}\left[\frac{C}{M_f-M_0}\right]$$ (37)

Moment calculation

Having defined the distribution function, the calculations of the regular and generalized moments are made according to Equations (1) and (2). For the regular moments, there is the following analytic solution:

$${\mu }_k=\overline{C}\left[\frac{\Gamma \left(\overline{A}\right)}{\Gamma \left(A\right)}\right]\left[\frac{1}{{\overline{B}}^{-k}}\right]\frac{\gamma \left(\overline{A},\overline{C}/\overline{B}\right)}{\Gamma \left(\overline{A}\right)}$$ (38)

where

$$\overline{A}=A+k$$

and γ(a, x) is the incomplete gamma function, defined as:

$$\gamma \left(a,x\right)=\underset{0}{\overset{x}{{\displaystyle \int }}}{t}^{a-1}\text{exp}\left(-t\right)dt$$ (39)

For the definition of the generalized moments, the Jacobi orthogonal polynomial family (Sack and Donovan, 1972) was chosen. This family has two parameters, α and β, which, when properly chosen, generate other polynomial families known in the literature.

The implementation of these polynomials consisted in computing the recursion coefficients a _k , b _k and c _k , used in Equation (3). These coefficients were generated by the method described by Press et al. (1992) for the interval [1, 1] and then modified to represent the orthogonal polynomial family in the interval [0, C].

Adiabatic flash calculation

Once the distribution function was discretized, the resulting composition was used for the adiabatic flash calculations, whose methodology can be found in Henley et al. (2011).

The methodology results in finding the composition of both vapor (x _i ^V ) and liquid (x _i ^L ) phases, the flash temperature (T ^flash ) as well as the vaporized fraction (γ) of the flash, with the feed stream and flash conditions given in Table 2.

Table 2
Adiabatic flash conditions.

	Distribution 1	Distribution 2
Feed temperature, TF	500 K	625 K
Feed pressure, PF	2 bar	3 bar
Flash pressure, Pflash	1 bar	1 bar

Conditions for the numerical analysis

Table 3 lists the methods compared for generating the quadrature and the type of moments employed by them. The analysis of the efficiency and robustness of these algorithms was carried out by calculating the computational cost and the errors in the reconstructed moments.

Table 3
Methods analyzed.

Method	Type of moments
LQMDA (Sack and Donovan, 1972)	Generalized
Chebyshev Algorithm (Chebyshev, 1858)	Generalized
LQMDA (Sack and Donovan, 1972)	Regular
Chebyshev Algorithm (Chebyshev, 1858)	Regular
GWA (Golub and Welsch, 1969)	Regular
PDA (Gordon, 1968) - PDA1	Regular
PDA (Lage, 2007) - PDA2	Regular

The computational cost was measured by the function clock, inside the file time.h (in C language). This function was called before and after each routine that computes the quadrature and, therefore, the reported CPU clocks are just for the quadrature computation.

The error of the reconstructed moments was calculated as the mean square of the relative errors (MSRE) of the moments between the analytical moments, given by Equations (1) or (2) (µ_k,dist ), and the reconstructed moments, given by Equations (4) or (5) (µ_k,reconst ), according to:

$$MSRE=\sqrt{\frac{1}{2n}{\displaystyle \sum }_{k=0}^{2n-1}{\left(\frac{{\mu }_{k,dist}-{\mu }_{k,reconst}}{{\mu }_{k,dist}}\right)}^2}$$ (40)

Test 1: Efficiency of the methods. The parameter C was kept equal to 1 and the computational cost and the MSRE were computed for all cases, varying the number of quadrature points from 3 to 20. Ten runs of the computer program were used to evaluate the computational cost, giving its mean value and standard deviation.

Test 2: Choice of distribution variable. The value of the parameter C was varied from 0.5 to 2.0 and its influence on the values of the 20 first reconstructed moments was analyzed for methods using either regular or generalized moments. Besides, for each value of C, the maximum number of quadrature points that can be computed by the method (n _max ) and the associated MSRE were evaluated as a measure of the method robustness.

Test 3: Choice of orthogonal polynomial family. Aiming at analysing the influence of the choice of the orthogonal polynomial family used to compute the generalized moments, the values of the parameters of the Jacobi polynomial were varied from -0.5 to 2.0. Some of the chosen parameter values generate polynomial families known in the literature, as shown in Table 4. The values of the 20 first reconstructed moments were compared and the maximum number of quadrature points achieved by the method (n _max ) and the associated MSRE were also evaluated.

Table 4
Parameters of the Jacobi polynomial families.

Polynomial Family	α	β
Chebyshev (1st order)	-0.5	-0.5
Chebyshev (2nd order)	0.5	0.5
Legendre	0	0

Test 4: Adiabatic flash solution. For this test, only the LQMDA using generalized moments was employed with C = 1 using the Jacobi polynomials with α = β = 2. The adiabatic flash was computed using different values for the number of quadrature points and the accuracy of the following properties of the streams were analyzed.

• Bubble point temperature of the feed stream at flash pressure (T ^bub );

• Dew point temperature of the feed stream at flash pressure (T ^dew );

• Flash temperature at flash pressure (T ^flash );

• Vapor fraction for the adiabatic flash (γ^flash );

• Mean molar mass for the feed, vapor and liquid streams (M ^F , M ^V , M ^L ).

The mean molar mass of a stream G is given by:

$${\overline{M}}^G={\displaystyle \sum }_{i=1}^n{M}_i^G{x}_i^G,G=F,L,V$$ (41)

These properties were compared to those obtained by using a uniform discretization of the feed molar fraction distribution using 10 to 10,000 pseudocomponents.

The moments for the liquid and vapor streams were computed from the flash results. Then, for each stream, a new characterization was calculated by computing a new Gauss-Christoffel quadrature. This was used to obtain the reconstructed moments of the stream and the corresponding MSRE values.

RESULTS

Test 1 results

The analysis of the efficiency and robustness of several algorithms for computing the Gauss-Christoffel quadrature (see Table 3) is given by the results shown in Tables 5 and 6 for their computational cost and in Tables 7 and 8 for the corresponding MSRE values, considering the truncated gamma distributions given in Table 1 for the molar fraction.

Table 5
Computational cost for the generation of the Gauss-Christoffel quadrature for the distribution 1, in clocks.

n	LQMDA	Cheb	LQMDA	Cheb	GWA	PDA2	PDA1
	gen	gen	reg	reg	reg	reg	reg
3	135.3 ± 1.2	8.7 ± 0.5	6.0 ± 0.0	3.9 ± 0.4	4.5 ± 0.6	16.1 ± 0.6	6.7 ± 0.5
4	141.8 ± 6.1	10.7 ± 0.5	7.8 ± 0.5	5.9 ± 0.4	7.1 ± 0.6	22.7 ± 1.8	11.8 ± 1.1
5	143.9 ± 12.3	13.1 ± 2.0	10.0 ± 0.9	7.9 ± 0.8	9.4 ± 1.0	29.8 ± 3.0	11.3 ± 1.0
6	149.1 ± 0.9	16.4 ± 0.6	13.4 ± 0.6	11.5 ± 0.6	16.1 ± 0.4	37.1 ± 1.2	14.8 ± 0.5
7	151.7 ± 1.4	19.1 ± 0.4	15.9 ± 0.4	1.4 ± 0.6	17.5 ± 0.6	46.0 ± 1.2	NA
8	155.4 ± 2.9	23.6 ± 0.7	22.5 ± 0.8	17.9 ± 0.8	24.0 ± 4.0	59.3 ± 5.4	NA
9	157.9 ± 13.5	27.4 ± 3.4	27.9 ± 6.1	22.3 ± 2.3	30.2 ± 5.1	65.3 ± 8.3	NA
10	167.9 ± 4.9	34.8 ± 0.5	30.3 ± 3.9	32.1 ± 6.8	34.8 ± 0.5	86.1 ± 10.0	NA
11	174.0 ± 6.6	42.4 ± 5.6	35.9 ± 4.0	34.7 ± 1.0	42.9 ± 0.4	98.6 ± 5.6	NA
12	185.0 ± 6.9	51.4 ± 3.2	42.3 ± 0.5	39.6 ± 5.0	59.9 ± 6.5	112.8 ± 8.7	NA
13	188.0 ± 6.8	52.5 ± 2.1	NA	NA	NA	NA	NA
14	198.3 ± 4.2	60.1 ± 0.4	NA	NA	NA	NA	NA
15	211.0 ± 8.0	68.3 ± 0.7	NA	NA	NA	NA	NA
16	218.9 ± 7.3	79.9 ± 4.5	NA	NA	NA	NA	NA
17	228.3 ± 6.0	88.4 ± 2.8	NA	NA	NA	NA	NA
18	236.3 ± 6.8	91.6 ± 7.0	NA	NA	NA	NA	NA
19	249.0 ± 7.2	108.8 ± 1.3	NA	NA	NA	NA	NA
20	259.2 ± 4.7	119.5 ± 1.2	NA	NA	NA	NA	NA

Table 6
Computational cost for the generation of the Gauss-Christoffel quadrature for the distribution 2, in clocks.

n	LQMDA	Cheb	LQMDA	Cheb	GWA	PDA2	PDA1
	gen	gen	reg	reg	reg	reg	reg
3	135.4 ± 1.1	8.4 ± 0.5	6.0 ± 0.0	3.9 ± 0.3	5.0 ± 0.0	16.5 ± 1.1	7.0 ± 0.0
4	142.4 ± 1.1	10.6 ± 0.5	8.0 ± 0.0	6.0 ± 0.0	7.3 ± 0.5	22.3 ± 1.3	11.5 ± 0.7
5	145.0 ± 1.3	12.9 ± 0.3	10.1 ± 0.3	8.4 ± 0.5	9.9 ± 0.3	30.9 ± 1.1	11.5 ± 0.5
6	148.3 ± 0.7	16.3 ± 0.5	13.5 ± 0.5	11.5 ± 0.5	16.3 ± 0.5	37.8 ± 1.1	14.9 ± 0.3
7	152.3 ± 1.5	20.0 ± 0.5	16.9 ± 0.3	17.2 ± 0.4	18.0 ± 0.5	46.2 ± 1.1	NA
8	156.7 ± 1.6	26.6 ± 10.0	22.3 ± 0.7	19.9 ± 5.3	22.6 ± 0.5	57.6 ± 0.7	NA
9	160.5 ± 1.7	27.9 ± 0.6	26.8 ± 0.6	23.4 ± 3.1	33.4 ± 5.5	69.7 ± 3.5	NA
10	167.0 ± 3.1	35.9 ± 0.3	29.8 ± 0.4	29.6 ± 0.5	35.1 ± 0.3	83.1 ± 0.7	NA
11	172.2 ± 4.1	42.5 ± 6.5	36.7 ± 7.5	35.9 ± 3.4	42.1 ± 0.3	99.7 ± 6.1	NA
12	179.2 ± 3.0	49.4 ± 5.9	NA	NA	NA	NA	NA
13	201.7 ± 6.4	57.2 ± 1.1	NA	NA	NA	NA	NA
14	200.4 ± 8.1	58.9 ± 3.9	NA	NA	NA	NA	NA
15	215.0 ± 5.4	68.0 ± 2.3	NA	NA	NA	NA	NA
16	224.8 ± 9.7	78.7 ± 0.7	NA	NA	NA	NA	NA
17	227.0 ± 5.6	83.8 ± 0.8	NA	NA	NA	NA	NA
18	231.9 ± 1.4	94.9 ± 0.6	NA	NA	NA	NA	NA
19	250.6 ± 7.8	111.0 ± 3.8	NA	NA	NA	NA	NA
20	252.4 ± 2.0	116.2 ± 1.1	NA	NA	NA	NA	NA

Table 7
MSRE in the reconstructed moments for the distribution 1.

n	LQMDA	Cheb	LQMDA	Cheb	GWA	PDA2	PDA1
	gen	gen	reg	reg	reg	reg	reg
3	7.02×10-16	1.16×10-15	9.01×10-16	3.16×10-15	3.25×10-13	9.40×10-16	9.40×10-16
4	2.31×10-15	4.73×10-15	2.42×10-15	1.87×10-15	2.97×10-4	7.78×10-16	7.78×10-16
5	9.74×10-15	1.10×10-14	6.60×10-15	1.09×10-15	2.42×10-5	4.16×10-15	4.16×10-15
6	9.24×10-15	1.28×10-14	1.07×10-14	7.35×10-16	1.83×10-6	3.40×10-15	3.40×10-15
7	2.61×10-14	1.04×10-14	1.07×10-15	5.70×10-15	1.33×10-7	3.81×10-15	NA
8	2.10×10-14	1.01×10-14	6.30×10-15	6.58×10-15	9.40×10-9	1.28×10-14	NA
9	5.70×10-14	7.03×10-14	4.21×10-15	3.87×10-15	6.50×10-10	3.90×10-16	NA
10	4.50×10-14	4.31×10-14	2.12×10-14	7.98×10-15	4.46×10-11	2.83×10-15	NA
11	1.80×10-13	1.65×10-13	2.57×10-15	8.30×10-15	2.76×10-12	1.76×10-14	NA
12	2.60×10-13	8.14×10-14	3.58×10-15	6.71×10-15	3.63×10-12	6.61×10-15	NA
13	5.77×10-13	2.59×10-13	NA	NA	NA	NA	NA
14	1.09×10-12	4.95×10-13	NA	NA	NA	NA	NA
15	2.38×10-13	1.69×10-13	NA	NA	NA	NA	NA
16	9.80×10-13	2.46×10-13	NA	NA	NA	NA	NA
17	4.78×10-13	5.53×10-13	NA	NA	NA	NA	NA
18	1.07×10-12	8.09×10-13	NA	NA	NA	NA	NA
19	1.77×10-13	1.69×10-13	NA	NA	NA	NA	NA
20	4.35×10-13	2.51×10-13	NA	NA	NA	NA	NA

Table 8
MSRE in the reconstructed moments for the distribution 2.

n	LQMDA	Cheb	LQMDA	Cheb	GWA	PDA2	PDA1
	gen	gen	reg	reg	reg	reg	reg
3	1.73×10-14	3.73×10-15	1.99×10-15	4.17×10-16	1.32×10-13	2.61×10-15	2.61×10-15
4	4.64×10-14	3.78×10-14	2.56×10-15	1.09×10-15	1.12×10-4	2.84×10-15	2.84×10-15
5	4.18×10-14	4.21×10-14	1.19×10-15	1.19×10-15	8.70×10-6	4.23×10-15	4.23×10-15
6	2.20×10-14	1.54×10-14	5.35×10-15	4.24×10-15	6.37×10-7	3.84×10-15	3.84×10-15
7	2.51×10-14	3.57×10-14	2.38×10-15	3.80×10-15	4.50×10-8	3.44×10-15	NA
8	4.90×10-14	4.16×10-14	3.84×10-15	1.83×10-15	3.10×10-9	1.22×10-15	NA
9	3.61×10-14	5.33×10-14	6.85×10-15	2.20×10-15	2.11×10-10	6.61×10-15	NA
10	4.83×10-14	7.99×10-14	2.82×10-15	1.02×10-14	1.38×10-11	6.18×10-15	NA
11	8.69×10-14	1.47×10-13	3.25×10-15	3.06×10-15	4.56×10-13	4.43×10-15	NA
12	7.25×10-14	6.69×10-14	NA	NA	NA	NA	NA
13	1.21×10-13	7.38×10-14	NA	NA	NA	NA	NA
14	1.66×10-13	4.61×10-14	NA	NA	NA	NA	NA
15	6.66×10-14	4.69×10-14	NA	NA	NA	NA	NA
16	6.87×10-14	6.04×10-14	NA	NA	NA	NA	NA
17	2.54×10-13	2.58×10-13	NA	NA	NA	NA	NA
18	2.79×10-13	4.22×10-13	NA	NA	NA	NA	NA
19	4.45×10-13	4.07×10-13	NA	NA	NA	NA	NA
20	3.61×10-13	4.88×10-14	NA	NA	NA	NA	NA

As expected, the computational cost increases with the number of quadrature points, The methods using regular moments are less expensive than those that employed the generalized moments. This is expected due to the additional cost of computing the recursion coefficients for the orthogonal polynomials. However, this cost increase for the Chebyshev method is not large.

The PDA2 was the most expensive method among those using the regular moments. This was expected due to the increase of operations related to the computation using the mantissa-exponent format. However, this made the method more robust than the original PDA (PDA1), allowing it to compute quadrature rules with a larger number of points.

In relation to the robustness of the methods, it can be observed that PDA1 is the least robust, being able to generate the Gauss-Christoffel quadrature rule for, at most, 6 quadrature points. The other methods using regular moments have similar behaviors among themselves, being able to compute the quadrature rule up to 12 points for distribution 1 and 11 points for distribution 2.

When analyzing the MSRE, although there are some oscillations, it can be noted that the MSRE tends to increase with the number of quadrature points due to error accumulation. The only exception was the GWA, for which the MSRE decreased by orders of magnitude as n increases. The reason for this behavior is that this method requires the usage of an additional moment of order 2n, which is the main factor responsible for the MSRE values for this method, because the n−point quadrature can only exactly compute the first 2n moments. For this method, the increase of the number of quadrature points increased the accuracy of this extra moment.

The methods using generalized moments were able to obtain the Gauss-Christoffel quadrature rule for more than 20 points for both distributions, showing that these methods are more robust than those using regular moments. The MSRE for these methods were slightly higher than for the other methods. However, these errors are of the order of 10⁻¹³, which is still very small for generating the Gauss-Christoffel quadrature rule for a large number of quadrature points.

Similarly to the results of John and Thein (2012), the LQMDA and the Chebyshev method were found to be equivalent in robusteness and computational cost when both used the regular moment set. However, for the generalized moment set, the Chebyshev method is much faster. For instance, considering the largest n value (11) for which the LQMDA and the Chebyshev method were able to compute the quadrature using both the generalized and regular moment sets for both distributions, the increase in the computational cost related to the usage of the generalized moments is about 376% for LQMDA and only 20% for the Chebyshev method.

Test 2 results

The effect of the C value used in the definition of the distribution variable on the moments used to compute the quadrature rule is shown in Tables 9 and 10 for molar fraction distributions 1 and 2, respectively.

Table 9
First 20 generalized and regular moments for the distribution 1 for several C values in the distribution variable definition.

C	0.5		1.0		1.333		1.5		2.0
k	mk(P)	mk	mk(P)	mk	mk(P)	mk	mk(P)	mk	mk(P)	mk
0	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100
1	-1.35×100	1.38×10-1	-1.35×100	2.75×10-1	-1.35×100	3.67×10-1	-1.35×100	4.13×10-1	-1.35×100	5.51×10-1
2	1.34×100	2.73×10-2	1.34×100	1.09×10-1	1.34×100	1.94×10-1	1.34×100	2.45×10-1	1.34×100	4.36×10-1
3	-1.02×100	6.85×10-13	-1.02×100	5.48×10-2	-1.02×100	1.30×10-1	-1.02×100	1.85×10-1	-1.02×100	4.38×10-1
4	8.39×10-1	2.02×10-13	8.39×10-1	3.24×10-2	8.39×10-1	1.02×10-1	8.39×10-1	1.64×10-1	8.39×10-1	5.18×10-1
5	-5.85×10-1	6.71×10-4	-5.85×10-1	2.15×10-2	-5.85×10-1	9.03×10-2	-5.85×10-1	1.63×10-1	-5.85×10-1	6.87×10-1
6	5.12×10-1	2.41×10-4	5.12×10-1	1.54×10-2	5.12×10-1	8.65×10-2	5.12×10-1	1.76×10-1	5.12×10-1	9.87×10-1
7	-3.50×10-1	9.19×10-5	-3.50×10-1	1.18×10-2	-3.50×10-1	8.80×10-2	-3.50×10-1	2.01×10-1	-3.50×10-1	1.51×100
8	3.43×10-1	3.66×10-5	3.43×10-1	9.38×10-13	3.43×10-1	9.35×10-2	3.43×10-1	2.40×10-1	3.43×10-1	2.40×100
9	-2.25×10-1	1.51×10-5	-2.25×10-1	7.73×10-13	-2.25×10-1	1.03×10-1	-2.25×10-1	2.97×10-1	-2.25×10-1	3.96×100
10	2.50×10-1	6.38×10-6	2.50×10-1	6.53×10-13	2.50×10-1	1.16×10-1	2.50×10-1	3.77×10-1	2.50×10-1	6.69×100
11	-1.52×10-1	2.75×10-6	-1.52×10-1	5.64×10-13	-1.52×10-1	1.33×10-1	-1.52×10-1	4.88×10-1	-1.52×10-1	1.16×101
12	1.93×10-1	1.21×10-6	1.93×10-1	4.95×10-13	1.93×10-1	1.56×10-1	1.93×10-1	6.42×10-1	1.93×10-1	2.03×101
13	-1.05×10-1	5.37×10-7	-1.05×10-1	4.40×10-13	-1.05×10-1	1.85×10-1	-1.05×10-1	8.57×10-1	-1.05×10-1	3.61×101
14	1.55×10-1	2.42×10-7	1.55×10-1	3.96×10-13	1.55×10-1	2.21×10-1	1.55×10-1	1.16×100	1.55×10-1	6.49×101
15	-7.41×10-2	1.10×10-7	-7.41×10-2	3.59×10-13	-7.41×10-2	2.68×10-1	-7.41×10-2	1.57×100	-7.41×10-2	1.18×102
16	1.30×10-1	5.02×10-8	1.30×10-1	3.29×10-13	1.30×10-1	3.27×10-1	1.30×10-1	2.16×100	1.30×10-1	2.15×102
17	-5.20×10-2	2.31×10-8	-5.20×10-2	3.03×10-13	-5.20×10-2	4.01×10-1	-5.20×10-2	2.98×100	-5.20×10-2	3.97×102
18	1.11×10-1	1.07×10-8	1.11×10-1	2.81×10-13	1.11×10-1	4.95×10-1	1.11×10-1	4.15×100	1.11×10-1	7.35×102
19	-3.56×10-2	4.98×10-9	-3.56×10-2	2.61×10-13	-3.56×10-2	6.15×10-1	-3.56×10-2	5.79×100	-3.56×10-2	1.37×103

Table 10
First 20 generalized and regular moments for the distribution 2 for several C values in the distribution variable definition.

C	0.5		1.0		1.333		1.5		2.0
k	mk(P)	mk	mk(P)	mk	mk(P)	k	mk(P)	mk	mk(P)	mk
0	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100
1	-6.46×10-1	1.96×10-1	-6.46×10-1	3.92×10-1	-6.46×10-1	5.23×10-1	-6.46×10-1	5.89×10-1	-6.46×10-1	7.85×10-1
2	2.93×10-1	4.71×10-2	2.93×10-1	1.89×10-1	2.93×10-1	3.35×10-1	2.93×10-1	4.24×10-1	2.93×10-1	7.54×10-1
3	1.09×10-1	1.32×10-2	1.09×10-1	1.05×10-1	1.09×10-1	2.50×10-1	1.09×10-1	3.56×10-1	1.09×10-1	8.44×10-1
4	-2.17×10-2	4.14×10-13	-2.17×10-2	6.62×10-2	-2.17×10-2	2.09×10-1	-2.17×10-2	3.35×10-1	-2.17×10-2	1.06×100
5	1.25×10-1	1.42×10-13	1.25×10-1	4.53×10-2	1.25×10-1	1.91×10-1	1.25×10-1	3.44×10-1	1.25×10-1	1.45×100
6	4.10×10-2	5.18×10-4	4.10×10-2	3.31×10-2	4.10×10-2	1.86×10-1	4.10×10-2	3.77×10-1	4.10×10-2	2.12×100
7	7.69×10-2	1.99×10-4	7.69×10-2	2.55×10-2	7.69×10-2	1.91×10-1	7.69×10-2	4.35×10-1	7.69×10-2	3.26×100
8	6.49×10-2	7.96×10-5	6.49×10-2	2.04×10-2	6.49×10-2	2.03×10-1	6.49×10-2	5.22×10-1	6.49×10-2	5.22×100
9	6.89×10-2	3.28×10-5	6.89×10-2	1.68×10-2	6.89×10-2	2.23×10-1	6.89×10-2	6.46×10-1	6.89×10-2	8.61×100
10	6.83×10-2	1.39×10-5	6.83×10-2	1.42×10-2	6.83×10-2	2.52×10-1	6.83×10-2	8.20×10-1	6.83×10-2	1.46×101
11	6.89×10-2	5.98×10-6	6.89×10-2	1.23×10-2	6.89×10-2	2.89×10-1	6.89×10-2	1.06×100	6.89×10-2	2.51×101
12	6.91×10-2	2.62×10-6	6.91×10-2	1.07×10-2	6.91×10-2	3.38×10-1	6.91×10-2	1.39×100	6.91×10-2	4.40×101
13	6.94×10-2	1.16×10-6	6.94×10-2	9.53×10-13	6.94×10-2	4.00×10-1	6.94×10-2	1.85×100	6.94×10-2	7.81×101
14	6.97×10-2	5.22×10-7	6.97×10-2	8.55×10-13	6.97×10-2	4.78×10-1	6.97×10-2	2.50×100	6.97×10-2	1.40×102
15	7.00×10-2	2.36×10-7	7.00×10-2	7.75×10-13	7.00×10-2	5.78×10-1	7.00×10-2	3.39×100	7.00×10-2	2.54×102
16	7.02×10-2	1.08×10-7	7.02×10-2	7.08×10-13	7.02×10-2	7.03×10-1	7.02×10-2	4.65×100	7.02×10-2	4.64×102
17	7.05×10-2	4.96×10-8	7.05×10-2	6.51×10-13	7.05×10-2	8.62×10-1	7.05×10-2	6.41×100	7.05×10-2	8.53×102
18	7.07×10-2	2.30×10-8	7.07×10-2	6.02×10-13	7.07×10-2	1.06×100	7.07×10-2	8.90×100	7.07×10-2	1.58×103
19	7.09×10-2	1.07×10-8	7.09×10-2	5.60×10-13	7.09×10-2	1.32×100	7.09×10-2	1.24×101	7.09×10-2	2.93×103

As can be seen, the range of values for the 20 first regular moments largely varies with the value of C. For instance, for C = 0.5, this range span 8-9 orders of magnitude and for C = 4/3 it spans just 2-3 orders of magnitude. This corroborates the value of 4/3 suggested by Lage (2007) to mitigate the numerical truncation errors in the quadrature rule computation, which involves the subtraction of products of these moments. Considering just the analyzed values, C = 4/3 and 1.5 are the best choice for distributions 1 and 2, respectively. Therefore, although it is not worthwhile to seek for a best C value, as it depends on the distribution, a choice of C within [4/3, 3/2] seems to be a good rule of thumb.

On the other hand, the value of C had no effect in the range of values of the generalized moments. This was expected because the change of variable made in the distribution function was also carried out in the orthogonal polynomials to maintain their orthogonality in the desired interval.

The effect of the distribution variable definition on the robustness and the efficiency of the methods can be analyzed from the maximum number of quadrature points for which the methods were able to generate the quadrature rule and the corresponding MSRE values, which are shown in Tables 11 and 12 for molar fraction distributions 1 and 2, respectively.

Table 11
Maximum number of quadrature points obtained for the Gauss-Christoffel quadrature for the distribution 1 for several C values in the distribution variable definition and the corresponding MSRE values.

C	LQMDA-gen		Cheb-gen		LQMDA-reg		Cheb-reg		GWA		PDA2		PDA1
	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE
0.500	85	5.57×10-11	85	2.42×10-11	12	3.58×10-15	12	6.71×10-15	12	3.63×10-12	12	6.61×10-15	6	3.40×10-15
0.700	85	2.29×10-11	85	3.14×10-12	11	7.57×10-15	11	9.90×10-16	11	5.21×10-15	11	4.42×10-12	6	5.13×10-15
0.900	85	3.43×10-12	85	4.32×10-11	11	1.07×10-14	11	8.01×10-15	11	4.92×10-12	11	3.36×10-15	6	6.64×10-15
1.000	85	5.57×10-11	85	2.42×10-11	12	3.58×10-15	12	6.71×10-15	12	3.63×10-12	12	6.61×10-15	6	3.40×10-15
1.100	85	3.35×10-11	85	6.93×10-12	11	1.91×10-14	11	1.02×10-14	11	4.45×10-12	11	1.20×10-15	6	2.17×10-15
1.200	85	1.95×10-11	85	2.90×10-12	11	3.74×10-15	11	9.52×10-15	11	4.07×10-12	11	1.49×10-14	6	3.73×10-15
1.300	85	2.89×10-11	85	3.13×10-11	12	4.87×10-15	12	2.79×10-15	12	6.61×10-12	12	5.05×10-15	6	3.74×10-15
1.333	85	3.42×10-12	85	1.83×10-11	11	9.82×10-15	11	1.19×10-14	11	4.78×10-12	11	3.85×10-15	7	4.25×10-12
1.500	85	2.60×10-11	85	2.96×10-12	11	7.79×10-15	11	5.38×10-15	11	4.53×10-12	11	1.69×10-14	7	8.43×10-15
1.750	85	1.93×10-11	85	4.96×10-12	11	9.17×10-15	11	9.79×10-15	11	4.13×10-12	11	1.78×10-14	7	1.20×10-14
2.000	85	5.57×10-11	85	2.42×10-11	12	3.58×10-15	12	6.71×10-15	12	3.63×10-12	12	6.61×10-15	7	3.81×10-15
2.500	85	6.55×10-12	85	9.60×10-12	11	4.17×10-15	11	8.29×10-15	11	2.89×10-12	11	2.07×10-14	7	1.05×10-14
3.000	85	2.60×10-11	85	2.96×10-12	11	7.79×10-15	11	5.38×10-15	11	4.53×10-12	11	1.69×10-14	8	4.66×10-15

Table 12
Maximum number of quadrature points obtained for the Gauss-Christoffel quadrature for the distribution 2 for several C values in the distribution variable definition and the corresponding MSRE values.

C	LQMDA-gen		Cheb-gen		LQMDA-reg		Cheb-reg		GWA		PDA2		PDA1
	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE	nmax	MSRE
0.500	84	4.90×10-12	84	9.46×10-12	11	3.25×10-15	11	3.06×10-15	11	4.56×10-13	11	4.43×10-15	6	3.84×10-15
0.700	84	1.48×10-13	84	2.73×10-12	11	1.20×10-14	11	1.25×10-15	11	7.95×10-14	11	1.25×10-14	6	1.87×10-15
0.900	84	8.41×10-12	84	7.04×10-12	11	3.42×10-15	11	5.93×10-15	11	9.65×10-13	11	9.19×10-16	6	1.09×10-15
1.000	84	4.90×10-12	84	9.46×10-12	11	3.25×10-15	11	3.06×10-15	11	4.56×10-13	11	4.43×10-15	6	3.84×10-15
1.100	84	4.76×10-12	84	5.30×10-12	11	2.74×10-15	11	4.42×10-15	11	2.20×10-13	11	1.15×10-14	6	6.78×10-15
1.200	84	3.85×10-12	84	4.70×10-12	11	6.56×10-15	11	1.49×10-15	11	1.49×10-12	11	7.28×10-15	6	2.77×10-15
1.300	84	2.84×10-12	84	1.72×10-11	11	3.01×10-15	11	1.22×10-15	11	2.33×10-12	11	1.09×10-14	7	7.56×10-15
1.333	84	8.95×10-12	84	7.35×10-12	11	7.17×10-15	11	1.37×10-14	11	1.44×10-12	11	3.50×10-15	7	1.56×10-15
1.500	84	1.01×10-11	84	1.00×10-11	11	7.12×10-15	11	2.42×10-15	11	3.98×10-13	11	2.32×10-15	7	4.82×10-15
1.750	84	7.83×10-12	84	2.91×10-12	11	2.83×10-15	11	2.75×10-14	11	1.29×10-12	11	5.51×10-15	7	4.33×10-15
2.000	84	4.90×10-12	84	9.46×10-12	11	3.25×10-15	11	3.06×10-15	11	4.56×10-13	11	4.43×10-15	7	3.44×10-15
2.500	84	1.11×10-11	84	5.93×10-12	11	2.05×10-15	11	1.53×10-15	11	1.09×10-12	11	6.01×10-15	7	2.40×10-15
3.000	84	1.01×10-11	84	1.00×10-11	11	7.12×10-15	11	2.42×10-15	11	3.98×10-13	11	2.32×10-15	8	1.64×10-15

As the order of magnitude of the generalized moments were not affected by the value of C used in the definition of the distribution variable, the maximum number of quadrature points does not change with it for the methods using such moments. Moreover, the corresponding MSRE values are basically independent of the C value.

On the other hand, the methods using the regular moments were affected by the value of C, notably PDA1, whose n _max value increased with the C value. This was not expected because the magnitude of the regular moments also increased with the the C value.

The other methods were able to calculate quadrature rules with 11-12 points without following any specific pattern. The MSRE values were around 10⁻¹⁴ - 10⁻¹⁵ for the PDA2, LQMDA-reg and Cheb-reg. The MSRE for GWA were higher, of the order of magnitude of 10⁻¹³ to 10⁻¹², but this can be explained by the error accumulation caused by the additional moment of order 2n _max .

Test 3 results

The results regarding the choice of the orthogonal polynomial family to generate the generalized moments were obtained by varying the values of the parameters α and β of the Jacobi polynomial family, The results for the distributions 1 and 2 are shown in Tables 13 and 14, respectively. It can be observed that the values of the 20 first generalized moments span a range that is reduced when the values of α and β were increased. For the tested values, it was observed that the generalized moments of the Legendre polynomials had the largest order of magnitude range, that is, |µ₀ ^(P) / µ_k ^(P) | ~ 10⁵ and 10⁶, k = 1, …, 20, for the distributions 1 and 2, respectively.

Table 13
Generalized moments for the distribution 1 for some orthogonal polynomial families.

α	2	2	1	1	0	1	0.5	0	-0.5
β	2	1	2	1	1	0	0.5	0	-0.5
k	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)
0	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100
1	-1.35×100	-6.23×10-1	-1.62×100	-8.98×10-1	-1.17×100	-1.74×10-1	-6.74×10-1	-4.49×10-1	-2.25×10-1
2	1.34×100	3.33×10-1	1.68×100	5.05×10-1	7.86×10-1	-1.13×10-1	2.11×10-1	1.84×10-13	-1.24×10-1
3	-1.02×100	-4.60×10-2	-1.39×100	-1.79×10-1	-3.64×10-1	1.41×10-1	2.88×10-2	1.29×10-1	1.49×10-1
4	8.39×10-1	4.74×10-2	1.07×100	5.10×10-2	1.20×10-1	-5.86×10-2	-7.17×10-2	-9.51×10-2	-6.62×10-2
5	-5.85×10-1	4.38×10-2	-7.95×10-1	5.11×10-13	-2.25×10-2	2.85×10-2	5.58×10-2	4.22×10-2	8.73×10-13
6	5.12×10-1	2.76×10-2	6.12×10-1	-3.89×10-13	-4.78×10-13	3.33×10-4	-2.99×10-2	-1.30×10-2	8.39×10-13
7	-3.50×10-1	4.21×10-2	-4.70×10-1	1.18×10-2	8.56×10-13	4.67×10-13	1.73×10-2	2.34×10-13	-8.33×10-13
8	3.43×10-1	3.11×10-2	3.81×10-1	-1.91×10-13	-6.94×10-13	4.82×10-13	-8.28×10-13	3.54×10-4	4.82×10-13
9	-2.25×10-1	3.87×10-2	-3.05×10-1	7.30×10-13	4.97×10-13	3.06×10-13	5.97×10-13	-6.71×10-4	-2.51×10-13
10	2.50×10-1	3.31×10-2	2.58×10-1	2.15×10-4	-3.54×10-13	3.77×10-13	-2.87×10-13	5.15×10-4	1.21×10-13
11	-1.52×10-1	3.72×10-2	-2.13×10-1	4.95×10-13	2.57×10-13	2.79×10-13	2.81×10-13	-3.49×10-4	-6.91×10-4
12	1.93×10-1	3.41×10-2	1.86×10-1	1.02×10-13	-1.92×10-13	3.02×10-13	-1.23×10-13	2.35×10-4	3.69×10-4
13	-1.05×10-1	3.66×10-2	-1.56×10-1	3.70×10-13	1.47×10-13	2.50×10-13	1.61×10-13	-1.63×10-4	-2.59×10-4
14	1.55×10-1	3.47×10-2	1.40×10-1	1.31×10-13	-1.15×10-13	2.55×10-13	-5.76×10-4	1.15×10-4	1.47×10-4
15	-7.41×10-2	3.62×10-2	-1.19×10-1	2.97×10-13	9.19×10-4	2.23×10-13	1.04×10-13	-8.21×10-4	-1.21×10-4
16	1.30×10-1	3.50×10-2	1.09×10-1	1.38×10-13	-7.51×10-4	2.22×10-13	-2.81×10-4	5.89×10-4	6.72×10-4
17	-5.20×10-2	3.61×10-2	-9.36×10-2	2.51×10-13	6.32×10-4	2.02×10-13	7.36×10-4	-3.99×10-4	-6.26×10-4
18	1.11×10-1	3.52×10-2	8.71×10-2	1.34×10-13	-5.47×10-4	1.96×10-13	-1.44×10-4	2.63×10-4	3.14×10-4
19	-3.56×10-2	3.61×10-2	-7.51×10-2	2.22×10-13	4.96×10-4	1.84×10-13	5.71×10-4	-1.20×10-5	-3.27×10-4

Table 14
Generalized moments for the distribution 2 for some orthogonal polynomial families.

α	2	2	1	1	0	1	0.5	0	-0.5
β	2	1	2	1	1	0	0.5	0	-0.5
k	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)	mk(P)
0	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100	1.00×100
1	-6.46×10-1	-3.82×10-2	-1.04×100	-4.31×10-1	-8.23×10-1	1.77×10-1	-3.23×10-1	-2.15×10-1	-1.08×10-1
2	2.93×10-1	-1.03×10-1	5.43×10-1	-5.73×10-2	1.77×10-1	-2.54×10-1	-1.63×10-1	-2.23×10-1	-2.36×10-1
3	1.09×10-1	2.23×10-1	-6.99×10-2	1.95×10-1	1.50×10-1	9.30×10-2	1.89×10-1	1.62×10-1	1.21×10-1
4	-2.17×10-2	4.03×10-2	-6.92×10-2	-7.58×10-2	-1.43×10-1	5.18×10-2	-4.95×10-2	-1.25×10-2	2.20×10-2
5	1.25×10-1	6.92×10-2	9.09×10-2	2.40×10-2	5.19×10-2	-2.39×10-2	-1.32×10-2	-3.66×10-2	-4.61×10-2
6	4.10×10-2	8.79×10-2	-3.66×10-2	2.04×10-2	-3.21×10-4	2.37×10-2	2.41×10-2	2.38×10-2	1.78×10-2
7	7.69×10-2	6.75×10-2	2.64×10-2	-1.67×10-13	-1.12×10-2	9.28×10-13	-9.00×10-13	-6.10×10-13	6.73×10-4
8	6.49×10-2	7.74×10-2	5.01×10-4	1.28×10-2	7.96×10-13	6.28×10-13	4.42×10-13	-1.04×10-13	-4.48×10-13
9	6.89×10-2	7.32×10-2	8.27×10-13	4.97×10-13	-3.68×10-13	9.14×10-13	1.26×10-13	1.83×10-13	2.46×10-13
10	6.83×10-2	7.43×10-2	5.39×10-13	7.00×10-13	1.34×10-13	6.30×10-13	6.78×10-4	-1.05×10-13	-8.72×10-4
11	6.89×10-2	7.39×10-2	5.58×10-13	5.71×10-13	-4.08×10-4	6.59×10-13	1.34×10-13	4.26×10-4	5.05×10-5
12	6.91×10-2	7.39×10-2	5.07×10-13	5.50×10-13	1.09×10-4	5.82×10-13	7.11×10-4	-1.39×10-4	2.23×10-5
13	6.94×10-2	7.39×10-2	4.79×10-13	5.09×10-13	-2.53×10-5	5.48×10-13	8.93×10-4	3.92×10-5	-9.16×10-5
14	6.97×10-2	7.39×10-2	4.49×10-13	4.78×10-13	5.53×10-6	5.09×10-13	6.82×10-4	-9.34×10-6	4.95×10-6
15	7.00×10-2	7.40×10-2	4.25×10-13	4.50×10-13	-2.12×10-7	4.78×10-13	6.84×10-4	2.57×10-6	-3.92×10-5
16	7.02×10-2	7.40×10-2	4.02×10-13	4.25×10-13	6.74×10-7	4.50×10-13	5.93×10-4	2.45×10-7	-1.24×10-5
17	7.05×10-2	7.41×10-2	3.83×10-13	4.03×10-13	1.39×10-6	4.25×10-13	5.65×10-4	1.04×10-6	-2.07×10-5
18	7.07×10-2	7.41×10-2	3.65×10-13	3.83×10-13	1.44×10-6	4.03×10-13	5.12×10-4	1.41×10-6	-1.27×10-5
19	7.09×10-2	7.42×10-2	3.49×10-13	3.65×10-13	2.40×10-6	3.83×10-13	4.83×10-4	1.93×10-6	-1.36×10-5

The results for the n _max and the MSRE values are shown in Tables 15 and 16 for the distributions 1 and 2, respectively. The computations using the Chebyshev polynomials of 1^st order (α = β = −0.5) had the lowest value of n _max for distribution 2 and presented the highest order of magnitude for the MSRE for both distributions. The Legendre polynomials (α = β = 0) also showed large values for the MSRE for both distributions. The Jacobi polynomials with α = β = 1 or α = β = 2 showed low MSRE values, making them good choices for the cases analyzed.

Table 15
Maximum number of quadrature points obtained for the Gauss-Christoffel quadrature for the distribution 1 for some orthogonal polynomial families and the corresponding MSRE values.

Type		LQMDA		Cheb
		Gen		gen
α	β	nmax	MSRE	nmax	MSRE
2.0	2.0	85	5.57×10-11	85	2.42×10-11
2.0	1.0	85	4.31×10-12	85	3.25×10-12
1.0	2.0	85	2.86×10-11	85	1.70×10-11
1.0	1.0	85	9.07×10-12	85	5.47×10-11
0.0	1.0	85	4.13×10-11	85	7.26×10-12
1.0	0.0	85	5.46×10-12	85	1.44×10-11
0.5	0.5	85	2.99×10-11	85	1.52×10-11
0.0	0.0	85	2.36×10-11	85	3.90×10-11
-0.5	-0.5	85	2.16×10-10	85	4.82×10-10

Table 16
Maximum number of quadrature points obtained for the Gauss-Christoffel quadrature for the distribution 2 for some orthogonal polynomial families and the corresponding MSRE values.

Type		LQMDA		Cheb
		Gen		gen
α	β	nmax	MSRE	nmax	MSRE
2.0	2.0	84	4.90×10-12	84	9.46×10-12
2.0	1.0	84	2.32×10-12	84	9.62×10-13
1.0	2.0	84	1.02×10-12	84	1.31×10-11
1.0	1.0	84	3.61×10-12	84	5.97×10-12
0.0	1.0	84	1.51×10-10	84	3.77×10-10
1.0	0.0	83	3.60×10-10	84	2.74×10-10
0.5	0.5	84	1.56×10-11	84	3.31×10-12
0.0	0.0	84	6.85×10-10	84	3.51×10-10
-0.5	-0.5	81	1.72×10-10	82	1.15×10-10

Test 4 results

For the adiabatic flashes of both distributions, according to the conditions given in Table 2, the results for the MSRE values and for some properties of the equilibrium streams for several numbers of discretization points are shown in Tables 17 and 18 for the distributions 1 and 2, respectively. It should be noted that the MSRE values for the liquid and vapor streams are those computed after the re-characterization of each stream, being similar to the MSRE results shown in Section 4.1 for the feed distributions.

Table 17
Stream properties for the adiabatic flash for distribution 1.

n	Tbub	Tdew	Tflash	gflash	$\overline{M}$ F	$\overline{M}$ V	$\overline{M}$ L	MSREF	MSREL	MSREV
	Gauss-Christoffel Quadrature
3	462.992	539.071	478.222	0.32949	155.087	131.721	166.568	1.04×10-15	2.17×10-15	4.05×10-15
4	461.939	539.336	478.312	0.31792	155.087	132.229	165.741	1.86×10-15	5.28×10-15	7.28×10-15
5	461.805	539.342	478.122	0.31756	155.087	132.281	165.699	7.51×10-15	9.61×10-15	4.99×10-15
6	461.791	539.342	478.121	0.31780	155.087	132.247	165.726	3.03×10-14	8.34×10-14	1.84×10-15
7	461.790	539.342	478.129	0.31779	155.087	132.248	165.725	1.50×10-14	8.21×10-14	2.49×10-14
8	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	1.80×10-14	8.10×10-14	4.30×10-14
10	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	5.29×10-14	4.22×10-13	1.41×10-14
12	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	2.29×10-13	1.22×10-13	1.08×10-14
14	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	6.07×10-13	1.30×10-13	1.20×10-14
20	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	2.03×10-13	3.68×10-13	4.85×10-13
30	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	1.85×10-11	8.20×10-13	1.48×10-12
50	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	7.96×10-12	5.38×10-13	1.82×10-13
80	461.790	539.342	478.128	0.31779	155.087	132.250	165.724	2.02×10-11	2.16×10-12	1.65×10-12
	Uniform Discretization
10	470.467	543.272	476.038	0.14931	157.854	133.576	162.115	NA	NA	NA
20	463.704	540.803	477.821	0.29001	155.739	132.517	165.225	NA	NA	NA
40	462.225	539.971	478.068	0.31161	155.268	132.304	165.662	NA	NA	NA
80	461.895	539.635	478.116	0.31624	155.146	132.260	165.731	NA	NA	NA
100	461.857	539.573	478.121	0.31677	155.130	132.255	165.735	NA	NA	NA
200	461.809	539.454	478.127	0.31748	155.104	132.250	165.735	NA	NA	NA
400	461.796	539.398	478.128	0.31767	155.094	132.250	165.731	NA	NA	NA
800	461.793	539.370	478.128	0.31773	155.090	132.250	165.728	NA	NA	NA
1000	461.793	539.364	478.128	0.31774	155.090	132.250	165.727	NA	NA	NA
2000	461.792	539.353	478.128	0.31776	155.088	132.250	165.725	NA	NA	NA
4000	461.792	539.347	478.128	0.31776	155.088	132.250	165.725	NA	NA	NA
8000	461.791	539.345	478.128	0.31777	155.087	132.250	165.724	NA	NA	NA
10000	461.791	539.344	478.128	0.31777	155.087	132.250	165.724	NA	NA	NA

Table 18
Stream properties for the adiabatic flash for distribution 2.

n	Tbub	Tdew	Tflash	gflash	$\overline{M}$ F	$\overline{M}$ V	$\overline{M}$ L	MSREF	MSREL	MSREV
	Gauss-Christoffel Quadrature
3	555.174	654.723	595.723	0.42818	237.324	196.014	268.257	1.73×10-14	1.67×10-15	3.95×10-16
4	551.918	654.896	594.212	0.42635	237.324	197.314	267.060	4.64×10-14	4.38×10-15	4.35×10-15
5	551.263	654.899	594.138	0.42662	237.324	196.756	267.509	4.18×10-14	1.11×10-14	4.21×10-14
6	551.156	654.899	594.221	0.42623	237.324	196.843	267.396	2.20×10-14	2.21×10-14	2.43×10-14
7	551.141	654.899	594.203	0.42644	237.324	196.852	267.415	2.51×10-14	2.24×10-14	1.21×10-13
8	551.139	654.899	594.205	0.42641	237.324	196.844	267.417	4.90×10-14	1.70×10-14	5.45×10-13
10	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	4.83×10-14	5.08×10-14	5.64×10-13
12	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	7.25×10-14	2.05×10-13	1.20×10-12
14	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	1.66×10-13	3.30×10-13	1.96×10-12
20	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	3.61×10-13	1.81×10-13	1.79×10-12
30	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	4.68×10-13	6.02×10-14	2.82×10-12
50	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	1.76×10-12	4.23×10-13	1.34×10-11
80	551.139	654.899	594.205	0.42640	237.324	196.846	267.415	5.04×10-13	3.79×10-12	9.91×10-12
	Uniform Discretization
10	551.288	657.766	594.440	0.42681	237.954	196.675	268.691	NA	NA	NA
20	551.144	656.338	594.288	0.42502	237.718	196.798	267.966	NA	NA	NA
40	551.178	655.615	594.243	0.42559	237.527	196.826	267.683	NA	NA	NA
80	551.161	655.255	594.224	0.42599	237.426	196.836	267.548	NA	NA	NA
100	551.157	655.184	594.220	0.42607	237.405	196.838	267.522	NA	NA	NA
200	551.148	655.041	594.213	0.42624	237.365	196.842	267.468	NA	NA	NA
400	551.143	654.970	594.209	0.42632	237.344	196.844	267.442	NA	NA	NA
800	551.141	654.934	594.207	0.42636	237.334	196.845	267.428	NA	NA	NA
1000	551.140	654.927	594.207	0.42637	237.332	196.845	267.426	NA	NA	NA
2000	551.139	654.913	594.206	0.42639	237.328	196.845	267.420	NA	NA	NA
4000	551.139	654.906	594.206	0.42640	237.326	196.846	267.418	NA	NA	NA
8000	551.139	654.902	594.205	0.42640	237.325	196.846	267.416	NA	NA	NA
10000	551.139	654.901	594.205	0.42640	237.325	196.846	267.416	NA	NA	NA

It can be observed that the discretization using the Gauss-Christoffel quadrature was not only capable of representing well both mixture properties, but it also does this with a much smaller number of discretization points. For both distributions, 8 quadrature points were enough to accurately represent the properties of these mixtures. For uniform discretizations of these mixtures, Tables 17 and 18 show that the number of pseudocomponents has to be around 10⁴ in order to achieve similar accuracy in the computation of the mixture properties.

Above 8 quadrature points, the Gauss-Christoffel quadrature discretization obtained the same values for all the properties of the streams, showing the good convergence of this method. The usage of more than 8 quadrature points only increased the MSRE values as already discussed.

For the uniform discretization method, it can be seen that its convergence is very slow, which can be verified by the fact that the stream properties still vary for n > 1000. Due to the large number of pseudocomponents that is needed, their molar fraction became quite low. This leads to a large accumulation of truncation errors that precluded the computations for n > 10000 for both distributions.

CONCLUSION

Among the methods compared for the computing of the Gauss-Christoffel quadrature, we concluded that those using generalized moments are more robust. Although these methods require a higher computational effort than those using regular moments, the gain in robustness is worthwhile as the number of quadrature points can be larger than 80. Even though this is an advantage, it might not be necessary as 8 pseudocomponents were shown to be enough for representing the properties of the streams involved in the adiabatic flash for the two mixtures analyzed in this work. The only disadvantage in using generalized moments is their larger computational cost. The usage of the generalized moment set made the LQMDA almost five times slower, whereas the Chebyshev algorithm showed just a 20% increase in its computational cost. Therefore, the Chebyshev algorithm using generalized moments is recommended to be used in the QMoM.

Due to the small number of pseudocomponents needed for accurate results, the QMoM proved to be an efficient and computationally cheap method for performing thermodynamic calculations for continuous and multicomponent mixtures, as also pointed out by Lage (2007), Rodrigues et al. (2012) and Petitfrere et al. (2014).

ACKNOWLEDGMENTS

Paulo L. C. Lage acknowledges the financial support from CNPq, grants #456905/2014-6 and 305265/2015-6. Argimiro R. Secchi acknowledges the financial support from CNPq, grant #302893/2013-0. This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001.

REFERENCES

• Briesen, H., Marquardt, W. An adaptative multigrid method for steady-state simulation of petroleum mixtures separation processes. Industrial & Engineering Chemical Research, 42, 2334-2348, 2003. https://doi.org/10.1021/ie0206150
  » https://doi.org/10.1021/ie0206150
• Chebyshev, P. L. Sur les fractions continues. Journal de Mathématiques Pures et Appliquées, 3, 289-323, 1858.
• Cotterman, R. L., Prausnitsz, J. M. Flash calculations for continuous or semicontinuous mixtures using an equation of state. Industrial and Engineering Chemistry Process Design and Development, 24, 434-443, 1985. https://doi.org/10.1021/i200029a038
  » https://doi.org/10.1021/i200029a038
• Gautschi, W. Algorithm 726: ORTHPOL - a package of routines for generating orthogonal polynomials and gauss-type quadrature rules. ACM Transactions on Mathematical Software, 20, 21-62, 1994. https://doi.org/10.1145/174603.174605
  » https://doi.org/10.1145/174603.174605
• Golub, G. H., Welsch, J. H. Calculation of gauss quadrature rules. Mathematics of Computation, 23, 221-230, 1969. https://doi.org/10.2307/2004418
  » https://doi.org/10.2307/2004418
• Gordon, R. Error bounds in equilibrium statistical mechanics. Journal of Mathematical Physics, 9, 655-663, 1968. https://doi.org/10.1063/1.1664624
  » https://doi.org/10.1063/1.1664624
• Haynes, H. W., Matthews, M. A. Continuous-mixture vapor-liquid equilibria computations based on true boiling point distillations. Industrial & Engineering Chemical Research , 30, 1911-1915, 1991. https://doi.org/10.1021/ie00056a036
  » https://doi.org/10.1021/ie00056a036
• Henley, E. J., Seader, J. D., Roper, D. K. Separation process principles. Wiley, 2011. ISBN 9780470646113. https://Books.Google.Com.Br/Books?Id=Cuo8bwaacaaj
  » https://Books.Google.Com.Br/Books?Id=Cuo8bwaacaaj
• Huang, S. H., Radosz, M. Phase behavior of reservouir fluids III: molecular lumping and characterization. Fluid Phase Equilibria, 66, 1-21, 1991. https://doi.org/10.1016/0378-3812(91)85044-U
  » https://doi.org/10.1016/0378-3812(91)85044-U
• John, V., Thein, F. On the efficiency and robustness of the core routine of the quadrature method of moments (QMOM). Chemical Engineering Science, 75, 327-333, 2012. https://doi.org/10.1016/j.ces.2012.03.024
  » https://doi.org/10.1016/j.ces.2012.03.024
• Lage, P. L. C. The quadrature method of moments for continuous thermodynamics. Computers & Chemical Engineering, 31, 782-799, 2007. https://doi.org/10.1016/j.compchemeng.2006.08.005
  » https://doi.org/10.1016/j.compchemeng.2006.08.005
• Liu, J. L., Wong, D. S. H. Rigorous implementation of continuous thermodynamics using orthonormal polynomials. Fluid Phase Equilibra, 129, 113-127, 1997. https://doi.org/10.1016/S0378-3812(94)02678-5
  » https://doi.org/10.1016/S0378-3812(94)02678-5
• Marano, J., Holder, G. A general equation for correlating the thermophysical properties of n-parafins, n-olefins, and other homologous series. 3. asymptotic behavior correlations for thermal and transport properties. Industrial Engineering Chemistry Research, 36, 2399-2408, 1997. https://doi.org/10.1021/ie9605138
  » https://doi.org/10.1021/ie9605138
• Mcgraw, R. Description of aerosol dynamics by the quadrature method of moments. Aerosol Science and Technology, 27, 255-265, 1997. https://doi.org/10.1080/02786829708965471
  » https://doi.org/10.1080/02786829708965471
• Petitfrere, M., Nichita, D. V., Montel, F. Multiphase equilibrium calculations using the semicontinuous thermodynamics of hydrocarbon mistures. Fluid Phase Equilibria , 362, 365-378, 2014. https://doi.org/10.1016/j.fluid.2013.10.056
  » https://doi.org/10.1016/j.fluid.2013.10.056
• Press, W. H., Flannery, B. P., Teukolsky, S. A., Vetterling, W. T. Numerical recipes in C: The art of scientific computing. Cambridge University Press, 2 edition, 1992.
• Rodrigues, R. C., Ahon, V. R. R., Lage, P. L. C. An adaptative characterization scheme for the simulaion of multistage separation of continuous mixtures using the quadrature method of moments. Fluid Phase Equilibria , 318, 1-12, 2012. https://doi.org/10.1016/j.fluid.2012.01.015
  » https://doi.org/10.1016/j.fluid.2012.01.015
• Sack, R. A., Donovan, A. F. An algorithm for gaussian quadrature given modified moments. Numerische Mathematik, 18, 465-478, 1972. https://doi.org/10.1007/BF01406683
  » https://doi.org/10.1007/BF01406683
• Upadhyay, R. R. Evaluation of the use of the chebyshev algorithm with the quadrature method of moments for simulating aerosol dynamics. Journal of Aerosol Science, 44, 11-23, 2012. https://doi.org/10.1016/j.jaerosci.2011.09.005
  » https://doi.org/10.1016/j.jaerosci.2011.09.005
• Wheeler, J. C. Modified moments and gaussian quadratures. Journal of Mathematics, 4, 287-296, 1974. https://doi.org/10.1216/RMJ-1974-4-2-287
  » https://doi.org/10.1216/RMJ-1974-4-2-287
• Whitson, C., Brule, M. (Editors). Phase Behavior. SPE, Texas, 2000.

NOMENCLATURE

• A  Gamma distribution parameter
• A _n  Terminal matrix
• a _k  Recursion coefficient for an orthogonal polynomial family of order k
• B  Gamma distribution parameter
• b _k  Recursion coefficient for an orthogonal polynomial family of order k
• C  Upper limit for the change of variable
• c _k  Recursion coefficient for an orthogonal polynomial family of order k
• C _p  Heat capacity
• D  Diagonal matrix
• f  Molar fraction distribution function
• FN  Distribution normalization factor
• H  Enthalpy
• ∆H ^vap  Enthalpy of vaporization
• H _f,i  Enthalpy of formation
• I  Identity Matrix
• I  Distribution variable
• K  Equilibrium constant
• M  Molar mass
• M  Mean molar mass
• MSRE  Mean square error
• N  Number of quadrature points
• P  Pressure
• P_k (x)  Christoffel orthogonal polynomial of order k
• p_k (x)  Orthogonal polynomial of order k
• q _i  Eigenvector of the terminal matrix, associated with the eigenvalue x _i
• T  Temperature
• x _i  Abscissa of the Gauss-Christoffel quadrature

Greek letters
• α  Parameter for the Jacobi polynomials
• β  Parameter for the Jacobi polynomials
• η_k  Recursion coefficient of the Christoffel orthogonal polynomial of order k
• γ  Vaporized fraction
• Γ  Gamma function
• µ_k  Regular moment of order k
• µ_k ^(P)  Generalized moments of order k
• ω_i  Weight from the Gauss-Christoffel quadrature
• θ_k  Recursion coefficient of the Christoffel orthogonal polynomial of order k

Superscripts
• bub  Bubble point
• dew  Dew point
• F  Feed stream
• flash  Flash condition
• G  Feed, vapor or liquid stream
• gi  Ideal gas condition
• L  Liquid stream
• sat  Saturation condition
• V  Vapor stream

APPENDIX

Properties Correlations

In this appendix, the correlations found in the literature for the estimation of the properties of the generated pseudocomponents are listed.

Saturation Pressure (Huang and Radosz, 1991):

$$B_1=9.5046+0.016104M$$ (42)

$$B_2=\text{exp}\left(5.0237+0.72702ln\left(M\right)\right)$$ (43)

$$P_{sat}=100000\text{exp}\left(B_1-B_2/T\right)$$ (44)

where [P _sat ] = Pa, and [T] = K.

Ideal Gas Heat Capacity (Marano and Holder, (1997):

$$\begin{array}{ll}{C}_P^{gi}=\hfill & \left(-0.0919055+0.011308T\left(K\right)-{6.3792010}^{-6}{T}^2\left(K\right)+{1.4060510}^{-9}{T}^3\left(K\right)\right)\hfill \\ \hfill & \left(N_c+0.284370\right)R\hfill \end{array}$$ (45)

where N _c = (M - 2)/14, R is the ideal gas constant and [C _P ^gi ] = [R].

Ideal Gas Enthalpy of Formation (Marano and Holder, 1997):

$$h_f^{gi}=-8.3206\left(N_c+2.111890\right)R{T}_0$$ (46)

where [h _f ^gi ] = [R][T ₀] and T ₀ = 298.15 K.

Enthalpy of Vaporization (Marano and Holder, 1997):

$$\Delta H^{Vap}=\left(1+1.99516\left(N_c-0.112756\right)\right)R{T}_0$$ (47)

where [ΔH ^Vap ] = [R][T ₀].

Publication Dates

• Publication in this collection
  09 Dec 2019
• Date of issue
  Jul-Sep 2019

History

• Received
  03 Nov 2018
• Reviewed
  08 Mar 2019
• Accepted
  18 Apr 2019