An Adaptive Multigrid Method for Steady-State Simulation of

2334

Ind. Eng. Chem. Res. 2003, 42, 2334-2348

An Adaptive Multigrid Method for Steady-State Simulation of Petroleum Mixture Separation Processes Heiko Briesen* and Wolfgang Marquardt† Lehrstuhl fu¨ r Prozesstechnik, RWTH Aachen University, Turmstrasse 46, D-52056 Aachen, Germany

On the basis of improvements in the analytical techniques for petroleum fraction characterization, the available detail in compositional information will strongly increase in the future. The inclusion of this detailed characterization on a molecular level seems to be the only way to fundamentally enhance the predictive capabilities of models for petroleum processes. An algorithm is presented that allows the solution of such models without actually solving the model in terms of single-component balances. Instead, the model is formulated using a distribution function to represent the composition. Using an adaptive wavelet-Galerkin discretization in combination with a multigrid approach, the algorithm is capable of solving the model with high accuracy without actually solving the high-dimensional nonlinear model equations. By satisfaction of the needs for scalable model formulation and an error controlled solution at high compositional detail, the approach presented seems to provide a promising alternative to currently employed pseudocomponent methods. 1. Introduction Instead of considering a few (say roughly up to 10) chemical components in standard multicomponent mixtures, complex multicomponent mixtures may comprise several hundred components. From an economical point of view, the most important industry dealing with complex multicomponent mixtures is the petroleum industry. Although the technique presented in this work applies to a variety of multicomponent mixture processes, the focus of this paper is on the discussion and exemplification of the proposed technique to thermal hydrocarbon separation processes. In the case of petroleum mixtures, the term complex multicomponent mixtures is certainly justified. Crude oil may contain a huge number of distinct molecular species1 on the order of 106. This large number is a result of the different branchings of the carbon chains. Considering hydrocarbon molecules with 25 carbon atoms, the number of possible configurations already exceeds 600 million.2 To reduce the model dimensionality of such complex multicomponent mixture processes, pseudocomponent lumping is widely used in industry and academia.3,4 With lumping methods, the thermophysical behavior of the mixture is represented by means of a small set of fictitious pseudocomponents. However, like any other model reduction, the lumping method introduces errors compared to the original formulation. In practice, the lumping method and the number of pseudocomponents has to be chosen on the basis of experience rather than on sound physical or mathematical criteria, although it is well-known that the choice of pseudocomponents can have a significant influence on the quality of the phase equilibrium calculations3 to name just one example. Another shortcoming of established lumping methods is the fixed set of pseudocomponents. Such a representation is not capable of tracking a major change in * To whom correspondence should be addressed. E-mail: [email protected]. † E-mail: [email protected].

composition in time or space with sufficient accuracy. Because the representation is not adaptive, the number of pseudocomponents has to be chosen sufficiently large to satisfactorily represent any composition possibly arising throughout the process. This produces an unnecessary overhead in complexity because all pseudocomponents have to be specified, although some of them may only occur during a short period of time or in selected locations if at all. Despite its drawbacks, the pseudocomponent method is virtually exclusively used because of its quite straightforward usability in current simulation frameworks. All models for standard chemical mixtures can be readily used by treating the pseudocomponents just as actual chemical species. Another way of representing the mixture is the use of a continuous distribution function.5 However, the continuous model demands the formulation of a discretization scheme to facilitate numerical solution. In fact, the use of a collocation discretization with distinct components selected as collocation points leads to a problem formulation that is fully equivalent to the pseudocomponent approach.6 Because of this, pseudocomponent lumping can merely be seen as a special, very restricted type of discretization. The necessity of a formal discretization scheme hinders straightforward implementation in the current simulation frameworks. However, it increases flexibility and creates the possibility of remedying the above-mentioned problems inherently associated with pseudocomponent methods. The structure of the paper is outlined in the following. The continuous modeling approach will be discussed in some more detail in the following section focusing on the use of a distribution function for composition representation of a mixture. By means of a simple example model, it will be shown how modeling of thermal separation processes is, in principle, affected by this way of representation. An abstraction of this simple example will provide us with the general problem formulation for which the subsequently reported algorithms are applicable.

10.1021/ie0206150 CCC: $25.00 © 2003 American Chemical Society Published on Web 05/03/2003

Ind. Eng. Chem. Res., Vol. 42, No. 11, 2003 2335

Figure 1. Upper part: discrete (left) and continuous (right) representation of a composition. Lower part: discrete (left) and continuized (right) representation. The continuous distribution function does not have to be smooth.

The solution method to be presented is based on a wavelet-Galerkin discretization of the continuous problem formulation. In section 3, some of the basic properties of the Haar wavelet basis functions are given. Also the application of the wavelet-Galerkin discretization procedure to the arising type of equations will be discussed. The proposed multigrid algorithm is then presented in section 4. Besides the basic framework of the algorithm, the means for further tuning of the algorithm are also discussed in light of the presented example. By investigation of some exemplary results, the performance of the algorithm is discussed in section 5. The paper closes with a concluding summary also discussing future research perspectives. 2. Continuous Mixture Representation and Modeling 2.1. Continuous Distribution Function. Before the pseudocomponent method gained popularity in the late 1960s,7 the continuous description was the prevailing way of characterization.8 The primary idea of the continuous mixture representation was to assume that the number of chemical species present in a petroleum (or other complex multicomponent) mixture is so large that it can be considered as a continuous rather than a discrete distribution, as shown in the upper part of Figure 1. The value of the molar fraction of a particular component class can always be determined by integrating the distribution F(ξ) over the corresponding interval of the characterizing variable ∆ξi:

∫∆ξ F(ξ) dξ ) xi i

(1)

For the rest of this work, only functions defined in the interval [0, 1] will be considered. This can be done without restriction of generality because every function defined in a finite interval can be scaled to the above interval by a simple linear transformation. There is large flexibility in the choice of the characterizing variable. Reasonable choices are, for example, the

relative volatility8 or the boiling temperature. However, a discrete component index for identifiable components may also be used. In the mid-1980s, many researchers5,9 have focused their attention on the continuous representation concept. The result of this development was a complete thermodynamic framework nowadays known as continuous thermodynamics. During implementation of the concepts of continuous thermodynamics, several technical problems such as flash calculations10 and multicomponent distillation11 have been addressed. A drawback of the above examples is that they use statistical distribution functions (e.g., Γ distribution) for the approximation of the distribution functions, which restricts the flexibility of the method and prevented its use in industrial practice. However, the general idea of representing a mixture composition by means of a distribution function is more general and not restricted to the continuous thermodynamics framework. Generally, any mixture with a finite number of components can be represented by a distribution function, like that depicted in the lower part of Figure 1. We can just use piecewise constant functions, where each function value corresponds to the discrete molar fraction of a component. The characterizing variable ξ is then no longer a true property coordinate of the components but just a component index. For example, to obtain the molar fraction of the second component, the integral in eq 1 over the interval [1/16, 2/ ] needs to be evaluated. Because of its foundation in 16 the discrete representation, the term continuized representation of a discrete mixture will be used. So, the term “continuous” here only states that for each value of ξ there is a corresponding value F(ξ). Therefore, a continuous distribution in the sense of this work does not have to be smooth. Note that for the continuized distribution function all property calculations can be performed in the standard discrete manner. Whether the function is smooth or highly irregular largely depends on the available data for characterization of the petroleum mixture. In the future, the analytical information will most likely comprise detailed data on a molecular level. However, currently the characterization of petroleum mixtures is usually done by means of distillation curves. To show the flexibility of the proposed approach, it will be shown that the success of the algorithm does not depend on the type of characterization. Note that that there is not any restriction on the order of the components. However, the success of the suggested algorithm relies on some similarity in the behavior of neighboring components. Because phase equilibrium is the governing phenomenon in distillation processes, an ordering by increasing normal boiling points is chosen. 2.2. Simple Continuous Flash Model. Modeling distillation processes with the continuous representation is quite similar to the standard way of discrete modeling. A simple steady-state flash as depicted in Figure 2 can, therefore, be modeled according to the following equations:

Material balance: 0 ) PfeedFF(ξ) - VFV(ξ) - LFL(ξ)

(2)

2336

Ind. Eng. Chem. Res., Vol. 42, No. 11, 2003

f[F ˆ (ξ),xˆ ,p] ) 0

(7)

g[O[F ˆ (ξ)],xˆ ,p] ) 0

(8)

The equations to be considered for discretization are the function-valued functions fi in the vector f. The arguments of these functions are the vectors F ˆ (ξ) of the state functions (distribution functions), xˆ of the state variables, and p of the parameters: Figure 2. Schematic for a continuously formulated steady-state flash.

Phase equilibrium: 0 ) FV(ξ) - K[FV(ξ),FL(ξ),T,p] FL(ξ)

(3)

Splitting factor: 0 ) RPfeed - V

(4)

Closure condition: 0)

∫01FV(ξ) dξ - 1

(5)

Total material balance: 0 ) Pfeed - V - L

(6)

Equations 2 and 3 relate the phase compositions represented by the distribution functions [FL(ξ) and FV(ξ)] and are therefore themselves distributed over the characterizing variable ξ. Note that there no longer is a vector equation setting up the material balance and the phase equilibrium for each pseudocomponent. Instead, there is one equation to be fulfilled for each value of the characterizing variable ξ. In this continuous modeling context, the usually discrete phase equilibrium constant also becomes a continuous function K(ξ). The set of scalar-valued equations is essentially the same as that in the discrete case. The only obvious difference is the replacement of the summation of all discrete molar fractions in the closure condition with the corresponding integration (5) over the whole range of ξ. Any process model available in the discrete representation can be rewritten without major effort in the continuous formulation. Basically, the only difference of a discrete and a continuous modeling approach is the replacement of all composition vectors by continuous distribution functions now representing the composition of the mixture in a phase. In the continuous formulation, there is not yet a measure for the complexity of the model like the number of pseudocomponents determining the total number of equations. The question of complexity of the continuous formulation remains unanswered until a discretization scheme to facilitate a numerical solution is introduced. The model complexity can then be considered as the number of equations emerging after discretization of each of the equations (2) and (3). The problem of discretization of the continuously formulated models will be discussed in detail in section 3.2. 2.3. General Model Formulation. Though the generality of the employed method of weighted residuals does not imply this restriction, only the class of models which actually arises will be examined. Thus, a steadystate model like the one introduced in eqs 2-6 is considered. These classes of models can be written in the following general residual formulation:

f ) (f1, f2, ..., fn)T

(9)

ˆ 2, ..., F ˆ n) T F ˆ ) (F ˆ 1, F

(10)

xˆ ) (xˆ 1, xˆ 2, ..., xˆ m)T

(11)

p ) (p1, p2, ...)T

(12)

In eq 8 an operator O[‚] is applied to the vector F ˆ (ξ) of the state function. This operator O[‚] results in a vector of scalar values so that F ˆ (ξ) is not directly an argument of the vector function g. The proposed solution strategy will be developed for this general type of equation system. However, the knowledge of the system of interest will allow further tuning of the algorithm to enhance the performance of the method, as will be presented in section 4.2. 3. Multiscale Techniques One of the building blocks of the method proposed is the use of wavelets. Wavelets have found a widespread acceptance in many different areas of mathematics and engineering. In mathematics, wavelets have ripened to provide a general purpose tool for problems such as approximation, adaptive gridding, solution of integral equations, and many others.12,13 Besides the most popular application of wavelets for image compression, also process engineers start to exploit the benefits of wavelet techniques for the solution of process problems. The area of signal processing is the most common field where wavelet algorithms are used today.14,15 Making information more manageable and compressing the data without losing too much accuracy is the common objective of all of these applications. Wavelets seem to provide an excellent tool when data or signals show some inherent structure or correlation.16 In the case of multicomponent mixture process simulation, the compositional data to achieve a robust and quick solution of our problem need to be compressed. However, the level of accuracy in the calculation still has to be high. The prerequisite for a promising application of wavelets, namely, some coherence or correlation in the data, is certainly fulfilled. The use of pseudocomponents actually relies on the similarity of properties of similar components. Therefore, von Watzdorf et al.17 started to exploit wavelet-based technology for the simulation of multicomponent mixture thermal separation processes. 3.1. Basic Properties of the Haar Wavelets. Though the theoretical foundation for wavelet techniques was initiated in the late 1980s,18 function representation in a wavelet-analogous framework is much older. Haar basis functions have been known since the beginning of the last century19 though not under the name wavelets. Haar wavelets can be considered as the simplest kinds of wavelets. However, they still


show the remarkable property to localize and compress information. It heavily depends on the function to be represented which wavelet basis is most suited for the approximation. As introduced in section 2, the distributions to be represented can either be smooth functions from a continuous experiment like a true-boiling-point (TBP) analysis or highly irregular functions such as those from a gas chromatographic (GC) analysis. Besides Haar wavelets, von Watzdorf6 investigated the Hat wavelet basis for approximating the concentration distributions. Both bases are specializations of B-spline wavelets of general order on finite domains. von Watzdorf showed that Hat wavelets are only superior to Haar wavelets if the distribution to be approximated is smooth. However, in the long run the importance of highly irregular representations will increase in order to include molecular information in the modeling of complex multicomponent mixture processes.20,21 In this case, a nonsmooth representation provided by Haar wavelets seems to be the most adequate choice. Also, the discretization with Haar wavelets is simple to implement, and one can interpret the results similarly to pseudocomponents because of the representation being piecewise constant. This section introduces the notation for the used Haar wavelets. Restricting further statements to Haar wavelets may lack a certain mathematical rigor but provides more clarity for the intended application. A discussion of some basic properties that are particularly important for the presented approach (orthonormality and norm equivalence) is given in work by Briesen.22 For a more general introduction into wavelets, we refer to work by Dahmen.13 Let us consider an arbitrary function F ˆ (ξ) from L2 (the space of square integrable functions) defined on the domain 0 e ξ e 1, for which a suitable approximation is desired. The starting point for generating Haar wavelets is the scaling function:

φ(ξ) )

{

1, 0 e ξ e 1 0, else

}

(13)

To generate basis functions on different scales j, the refinement equation needs to be utilized:

φj,l ) 2 j/2φ(2 jξ - 1), with l ) 0, ..., 2j - 1 (14) For fixed j, this refinement equation will provide a set of basis functions covering the domain [0, 1]. For suitably high scales, any function F ˆ (ξ) can be represented by the function F(ξ), which is an expansion up to an arbitrary level of detail:

F ˆ (ξ) ≈ F(ξ) )

∑ cj,lφj,l(ξ)

(15)

(j,l)∈Γ

To ease the notation, the approximation F(ξ) with the full set of basis functions can also be written as a scalar product of a coefficient vector and a basis function vector:

F ˆ (ξ) ≈ F(ξ) ) cT‚Φ(ξ)

(17)

with

cT :) (..., cj,l, ...) Φ :) (..., φj,l, ...)T, ( j, l) ∈ Γ The mother wavelet, from which the Haar wavelet basis can be constructed, is given by

{

1, 0 e ξ < 0.5 ψ(ξ) ) -1, 0.5 e ξ < 1 0, else

}

(18)

Applying the refinement equation by analogy to the single-scale refinement equation (14) provides the hierarchical basis functions

ψj,k(ξ) ) 2j/2ψ(2jξ - k)

(19)

Using the scaling function φ0,0 and the introduced multiscale basis functions ψj,k, any function on the interval [0, 1] can be approximated up to an arbitrary level of detail depending on the value of the maximum scale ˆj: ˆj 2 j-1

F ˆ (ξ) ≈ F(ξ) ) c0,0φ0,0(ξ) +

∑ ∑ dj,kψj,k(ξ) j)0 k)0

(20)

To ease notation, we introduce the definitions

ψ-1,0(ξ) :) φ0,0(ξ), d-1,0 :) c-1,0

(21)

With these definitions, eq 20 can be rewritten as

F(ξ) )

∑

dj,kψj,k(ξ)

(22)

(j,k)∈Λ

with the index pair set Λ given by

Λ :) {(-1, 0)} ∪

(23)

{( j, l): j ) 0, ..., ˆj, l ) 0, 1, ..., 2j - 1}

(24)

Using vector notation yields the even more compact form

F(ξ) ) dT‚Ψ(ξ)

(25)

of the so-called multiscale representation with

The index pair set Γ depends on the highest scale ˆjs actually used for the approximation: ˆ

Γ :) {(j,l): j ) ˆjs, l ) 0, 1, ..., 2 js - 1}

(16)

For example, the function shown in the lower right plot of Figure 1 can be written on the scale ˆjs ) 4 using 16 basis functions. The representation in eq 15 is referred to as a single-scale representation and the coefficients cj,i, which determine the shape of the function, are called single-scale coefficients.

dT :) (..., dj,k, ...) Ψ : ) (..., ψj,k, ...)T, (j, k) ∈ Λ

(26)

Note that both the multiscale and single-scale representations are fully equivalent:

F ˆ (ξ) ≈ F(ξ) ) dT‚Ψ ) cT‚Φ

(27)

If ˆj or ˆjs ) ˆj + 1 is chosen suitably large, this multiscale or single-scale representation F(ξ) of the

2338


original function F ˆ (ξ) can be considered as the best possible representation. This representation can be used as a reference for reduced approximations. Fast Wavelet Transformation (FWT). To change from a single-scale representation to a multiscale representation, the following simple linear transformation can be employed:

d ) Tc

∑

dj,kψj,k

(29)

(j,k)∈Λr

The approximation error can be measured by the norm of the difference of the full and reduced representations. Additionally, exploiting the norm identity property for the Haar wavelets yields

) ||F(ξ) - F r(ξ)||L2 )

Fi(ξ) ) diT‚Ψ

(28)

with c and d both being vectors of the dimension n ) ˆ 2js and T being a n × n matrix. However, the transformation is not actually implemented following this matrix multiplication. The hierarchical structure of the wavelets allows one to set up a much more efficient transformation algorithm. While the general linear transformation in eq 28 requires O(n2) operations, the FWT algorithm can be realized at the cost of O(n). Because the same is true for the backtransformation, switching to and from the wavelet representation is very efficient. For the details of the FWT algorithm, we refer to, e.g., work by Sweldens and Schro¨der.16 Function Approximation. If not all of the basis functions up to a certain maximum scale in Λ, but only a subset in Λr ⊂ Λ, are used, a reduced approximation of the distribution function is obtained:

F(ξ) ≈ Fr(ξ) )

constructed in section 3.1 are used. The basis functions of these spaces are now used as approximating and projection functions of the wavelet-Galerkin scheme. If the Haar wavelet basis is used, any of the functions Fi(ξ) from the vector of distribution functions F(ξ) can be represented according to eq 25 as

∑

dj,k2

(30)

(j,k)∈Λ\Λr

In the expression on the right, the squares of the multiscale coefficients of the neglected basis functions are summed up. It is reasonable to neglect the basis functions in Λ with the smallest corresponding coefficient in d to reduce the complexity with a minimum effect on accuracy. With the vector norm of the neglected coefficients, a direct measure for the quality of the approximation is available. 3.2. Wavelet-Galerkin Discretization. In the field of numerical mathematics, multiscale functions such as wavelets have been used to set up wavelet-Galerkin discretization schemes.23,24 For example, in chemical engineering the wavelet-Galerkin method has been applied to solve population balance equations.25,26 von Watzdorf et al.6,17,27 developed an error-controlled and adaptive wavelet-Galerkin method to solve continuously formulated multicomponent separation problems as they have been introduced in section 2. In contrast to pseudocomponents, their approach allows one to specify an error criterion instead of a fixed number of components. The algorithm then determines a discretization, which specifies the optimal detail of the composition representation to solve the problem with a prescribed accuracy at minimum computational cost. Although this work follows a completely different strategy for error control and adaptation (see section 4), the basic discretization, the wavelet-Galerkin method, is adopted. To set up the discretization scheme, the approximation spaces (multiscale and single scale) that have been

(31)

with the vector of multiscale distribution coefficients di and the vector of wavelet basis functions Ψ corresponding to the full index pair set of basis function Λ up to a certain maximum scale ˆj. To obtain the residual formulation, eq 31 can be introduced in one of the functions fi from the vector function f in eq 7 to obtain

ri(ξ) ) fi[(d1T‚Ψ, ..., dnT‚Ψ)T,x,p)

(32)

Because it is the objective to set up a Galerkin discretization, the same wavelet basis functions are also used as projection functions in the framework of the method of weighted residuals:

(

)

〈r1(ξ), Ψ〉 ) f˜(d,x,p) ) 0 l 〈rn(ξ), Ψ〉

(33)

This discretization has led to a nonlinear algebraic system of equations f˜ depending on the vector x of the state variables and the vector d of all of the distribution coefficients

d ) (d1T, d2T, ..., dnT)T

(34)

To close this set of algebraic equations (33), it is augmented by the scalar-valued equations g from eq 8 in which the operator O[‚] is evaluated using the current representation of the distribution functions:

g˜ (d,x,p) ) 0

(35)

The augmented system can now be solved using, for example, a Newton-type method. Having this system solved, the unknown distribution functions Fi(ξ) are also determined according to eq 31. If the maximum scale ˆj, which determines the set of basis functions according to Λ, is large enough, this representation is considered as a high-resolution solution. Note, in the case of a discrete characterization of the multicomponent mixture using analytical techniques, the maximum of ˆj compares to the total number of chemical components identified in the mixture. Similar to the approximation of functions, it is not the objective to determine the high-resolution solution exactly in terms of a full set of basis functions up to scale ˆj. As seen in the approximation example, the key to an efficient reduction of complexity is to select only a subset of basis functions for the approximation. In the discretization problem, this corresponds to the objective not to determine all distribution coefficients d but only the most relevant. Generally, the set of basis functions to be selected in order to obtain the best results is not known a priori. This is contrary to the approximation example, where the coefficients are already known. How to properly select the best basis functions is the question of adaptivity, which will be discussed in the following section.


For now it is assumed that a reduced set of basis functions has been obtained by some means

Ψr ) (..., ψj,k, ...)T, ( j, k) ∈ Λr

(36)

⊂ Λ, which would corresponding to the index pair set give a reasonably good approximation of the solution. Provided this set Λr is given, the wavelet-Galerkin discretization leads to Λr

( )

〈f1, Ψr 〉 ) f˜ r(dr,x,p) ) 0 l r 〈fn, Ψ 〉

(37)

For illustration purposes, the application of the discretization scheme to the equations of the flash example is given in the appendix. Note that the nonlinear system (37) has only as many equations as needed to determine the hopefully significantly smaller set of distribution coefficients

dr ) [(dr1)T, ..., (drn)T]T

(38)

T r r dri ) (..., di,( j,k), ...), ( j, k) ∈ Λ

(39)

with

From the solution of the nonlinear system (37) augmented with the algebraic equations r

r

g˜ (d ,x,p) ) 0

(40)

only reduced approximations F ri of the unknown distribution functions are obtained because only a reduced set of basis functions has been employed:

Fi(ξ) ) diT‚Ψ ≈ F ri (ξ) ) (dri )T‚Ψr

(41)

Single-Scale Discretizations. Alternatively, the system can be discretized using the single-scale basis functions. The discretization is straightforward, following the same path as that for the multiscale approximation. Using the full single-scale basis for function representation according to eq 17 and for projection in the discretization, we obtain

( )

〈fi, Φ〉 l ) f˜*(c,x,p) ) 0 〈fn, Φ〉

(42)

Besides the obvious difference of the discretized functions f˜ and f˜*, the only difference between the multiscale and single-scale cases is that the distribution coefficients, the actual unknowns of the discretized system, are not the multiscale coefficients d but the single-scale coefficients c. The vector of single-scale coefficients is set up by analogy to eq 34. However, because single-scale functions have been used for approximation and projection, the resulting discretization scheme is no longer a wavelet-Galerkin discretization. Nevertheless, there is no difference in the resulting function representations because both representations are fully equivalent. However, as it will be

discussed in section 4, there are cases where it is more suitable to use the single-scale functions for discretization. 4. Adaptive Multigrid Solution Algorithm Today, multigrid methods are among the most powerful tools known in numerical mathematics with a wide range of applications.28,29 The success of multigrid methods stems from the fact that the problem never has to be solved on the fine grid but only on relatively coarse grids with strongly reduced computational effort. Using an initial guess, only a correction to this guess is computed. This correction is calculated on a relatively coarse grid and is interpolated to the fine grid. By correcting the initial guess with the interpolated correction term, one obtains a new approximation of the solution. Performing iterative loops, the approximations approach the solution. Because this work is written from an application perspective, it neither makes any attempt to review the broad literature of multigrid methods nor gives any concern to the proof of convergence for the proposed scheme. In fact, a rigorous proof of convergence of the multigrid method for the problem class considered in this work is still lacking. 4.1. Basic Multigrid Algorithm. Before discussion of the algorithm, first some notations are set up. The exact solution F ˆ (ξ) for any of the state functions of the problem in eqs 7 and 8 can be written as

F ˆ (ξ) ) Fappr(ξ) + ∆F(ξ)

(43)

where Fappr(ξ) is any approximation of the distribution and ∆F(ξ) is an unknown correction to obtain the exact solution. An approximation Fi(ξ) on the full grid results in

‚Ψ + ∆diT‚Ψ (44) F ˆ i(ξ) ≈ Fi(ξ) ) diT‚Ψ ) dapprT i Note that the representation Fi(ξ) gives the highresolution approximation used as a reference. If the approximation complexity is to be reduced further, only a reduced set of basis functions is used for the correction can be kept term ∆dri . However, the initial guess dappr i at the high resolution because it is already known:

‚Ψ + (∆dri )T‚Ψr F i(ξ) ) diT‚Ψ ≈ dapprT i

(45)

For clarity the proposed method will be discussed in light of the intended application to multicomponent mixture processes only. In section 3.2, a Galerkin discretization scheme has been set up to derive algebraic systems of equations for the approximate solution of the continuously formulated problem set up by eqs 7 and 8. It was shown that either full sets of projection functions Ψ or Φ or the reduced set of multiscale functions Ψr resulting in a model of reduced dimension may be employed. The proposed multigrid solution strategy now combines the basic ideas of multigrid methods and the wavelet-Galerkin discretization. Algorithm 1: Basic Multigrid Algorithm. Step 0. Set iteration counter z ) 0. Choose a reasonable initial guess for all of the state functions in the vector F0 in terms of the distribution coefficients d0 ) (d0,1T, ..., d0,nT)T on the full grid and the state variable vector x0. Also, choose a reduced discretization grid defined by Λ r.

2340


Step 1. In the first step, a reduced model for the correction of the current iterate needs to be set up and solved. If we now use an intermediate approximation

Fint z (ξ)

( )(

int T r T (dz,1 ) ‚Ψ dz,1T‚Ψ + ∆dz,1 ‚Ψr ) l ) l int T r T dz,nT‚Ψ + ∆dz,n ) ‚Ψ ‚Ψr (dz,n

)

(46)

in the Galerkin discretization employing the same set of basis functions for projection as for the approximation of the correction term, a nonlinear system of equations r and xint for the unknowns ∆dz,i z is obtained:

0)

(

) (

r int 〈f(Fint f˜r(∆drz,xint z (ξ),xz ,p), Ψ 〉 z ,p) ) int int r int g{O[Fz (ξ)],xz ,p} g˜ (∆dz,xz ,p)

)

(47)

with the correction coefficients arranged in the vector r r T , ..., ∆dz,n ) ∆drz ) (∆dz,1

(48)

and the vector xint denoting the intermediate state z variables. Note that we no longer solve for the distribution coefficients representing the state functions Fz. In contrast to the expositions in section 3.2, the system is r r T ) ∆dz,i ‚Ψr to formulated in terms of corrections ∆Fz,i the current iterate. Step 2. Solve the nonlinear system (47) with a suitable nonlinear solver (e.g., a Newton-type method). Step 3. Check an error criterion for the intermediate approximation. If the accuracy constraint is fulfilled, terminate; else, continue. This step will be detailed in section 4.2.2. Step 4. To interpolate the reduced correction term to the full grid, a full grid linearization of eqs 33 and 35 is used. The linearization is done at the intermediate approximation according to

0)

(

int f˜(dint z ,xz ,p) int int g˜ (dz ,xz ,p)

)

+

()| (

( )

∂

˜f(‚) g˜ (‚)

∂

d x

T

dz+1 - dint z xz+1 - xint z

)

(49)

z,int

Solving this linear system of equations for the unknowns dz+l and xz+l provides the new iterates of the approximation for the next loop. Though the computational cost to solve this linear system should be much lower than that for the original full grid nonlinear system of equations, it still would be considerably large. However, physical insight allows one to simplify the solution of the system, as will be discussed in section 4.2. Step 5. Set z :) z + 1. Go back to step 1 with the new full grid approximation. The above introduced algorithm can also be seen as a Newton-Raphson scheme with inner and outer iteration loops. The connection to Newton’s method becomes obvious if algorithm 1 is represented in a different, more descriptive way like in Figure 3. Starting with a full grid approximation of the solution (step 0), the reduced problem is set up (step 1). This reduced nonlinear system is solved with a Newton-Raphson method (step 2). This part is the first iterative process, the inner loop.

Figure 3. Representation of the proposed basic multigrid algorithm in terms of an inner and outer Newton-type algorithm.

This loop is performed multiple times until the reduced system is actually converged. With this reduced grid solution for the correction term, an intermediate full grid approximation is constructed. After the termination check (step 3), the linearization of the full grid system at the intermediate approximation just gives an equation that actually constitutes one full grid Newton step from the intermediate solution to the next iterate (step 4). This outer Newton step is performed only once. Finally, we return to the construction of the reduced model for the new reduced approximation of the error with respect to the current full grid approximation of the solution. However, this should not be mistaken for a Newton multigrid method, in which the nonlinear system is usually solved by the standard sequence of linearized problems while only the underlying linear problems are tackled with an iterative multigrid solver. 4.2. Tuning of the Basic Algorithm. In the basic algorithm, it was already mentioned that physical reasoning can be used to increase the efficiency of the method. In the following, this potential will be exploited by modifying and extending several parts of the basic algorithm. 4.2.1. Efficient Interpolation of the Reduced Grid Correction. First the tuning of the calculation of the high-resolution correction term (step 4 in algorithm 1) should be considered. Because the whole algorithm is iterative, providing us only with ap-


proximations of the solution, the system (49) need not to be solved exactly. Rather, physically motivated simplifications and rearrangements of the actual equations can be introduced to significantly reduce the complexity of this step. The continuously formulated model equations for a flash contain the flow variables as well as temperature and pressure. As the composition representation Fz is improved by each iteration, the state variables xz become less dependent on the actual composition and are therefore considered to be constant in the linearized equation. With constant state variables, only the continuously formulated equations have to be processed because for the state variables x no estimates need to be determined. To employ further physically motivated simplifications, it is much more convenient to linearize the system in the discretized single-scale formulation (42). The linearization of the single-scale formulation under the assumption of constant state variables leads to

1)

0)

Pint feed

F cj,l

-

V Vintcj,l

-

L Lintcj,l

(51)

Also, the vapor-liquid equilibrium has to be linearized. However, the vapor-liquid-phase equilibrium constant can be regarded to show only a weak dependence on the concentration. This is especially true for weakly polar petroleum mixtures. For evaluating the phase equilibrium constants, the concentration information of the intermediate approximation can be used: V int L - kj,l cj,l, ( j, l) ∈ Γ 0 ) cj,l

(52)

int are Note that the phase equilibrium constants kj,l known. They can directly be calculated using the intermediate approximation of the solution. With this simplification, the vapor-liquid equilibrium is also already linear, which makes formal linearization unnecessary. However, first formally linearizing the equations and then introducing the above simplification lead to the same results. The number of equations of this linear problem is ˆ 2‚2 js. To reduce this number, the vapor distribution coefficients cV can be eliminated simply by introducing eq 52 into eq 51. This decreases the size of the system ˆ to 2 js equations. The solution of the linearized system would not satisfy the closure condition

T

s

(53)

(j,l)∈Γ

for either phase L or V. One way to overcome this problem is to normalize the solution of the linearized system in order to force the fulfillment of eq 53. However, after this normalization, the material balance would no longer be satisfied. To include the closure condition in the set of equations, some relaxation is necessary. A reasonable way to do this is to allow the phase equilibrium equation to be violated. For this purpose, the phase equilibrium conint are scaled with a factor q: stants kj,l V int L V int L - kj,l cj,l + ): cj,l - qkj,l cj,l 0 ) cj,l

(54)

This corresponds to approximately relaxing the system pressure from the known value pint to the unknown pressure p. For illustration, consider an ideal vaporliquid equilibrium V ) cj,l

Instead of solving for the correction term ∆c, here the correction equation ∆cz ) cz+l - cint z is already introduced, giving rise to a problem formulation in terms of the new approximation cz+l for the next iteration. For the flash example, the material balance equation does not have to be formally linearized. Provided the state variables V and L are assumed to be constant, the material balance equations are already linear. So, for each (j, l) ∈ Γ, this results in (for clarity the iteration counter z is omitted in eqs 51-56)

∫01(cL/V) ‚Φ dξ ) 2-jˆ /2 ∑ cj,lL/V

/ pj,l

p

L cj,l )

/ pint pj,l

pint

p p

p

cL ) int j,l

int L kj,l cj,l

/ is the pure-component vapor pressure of the where pj,l component associated with the single-scale indices (j, l). If eqs 55 and 54 are compared, the relaxation can be considered as a variable replacement q ) pint/p. Linearizing eq 54 and eliminating cV from eq 51 result in the desired linear system. Augmenting this linearized system with the closure condition for the liquid phase yields

F int Lint int L Pfeedcj,l + V intkj,l cj,l ) (V intkj,l + Lint)cj,l + int L,int cj,l )q, ∀ ( j, l) ∈ Γ (55) (V intkj,l ˆ

1 ) 2-js/2

∑ cj,lL

(56)

(j,l)∈Γ

Having this linear system solved for the liquid-phase coefficients cL and the pressure correction parameter q, one can evaluate eq 54 to determine cV. Therefore, our initial task, namely, to determine a new highresolution iterate for the distribution coefficients, is accomplished. By elimination of the vapor concentration variables, introduction of simplifications for the phase equilibrium constant, and relaxation of the phase equilibrium equation, it was possible to set up a problem in which the problem size was reduced by roughly 50%. The solution of this smaller problem obeys the closure condition and also satisfies the material balance. Thus, step 4 in algorithm 1 can be replaced by the following algorithm. Algorithm 2: Flash-Specific Replacement for the Basic Multigrid Algorithm. Step 4.1. Use the FWT to transform the intermediate approximations of the multiscale distribution coeffor both phases to the single-scale coefficients dL/V,int z ficients:

cL/V,int ) T-1dL/V,int z

(57)

Step 4.2. Solve the linear system (55) and (56) for cL and q. Step 4.3. Calculate cV from eq 54. Step 4.4. Update the multiscale distribution coefficients by retransformation of the single-scale coef-

2342


ficients:

can be easily calculated by evaluating the discretized residuals on the full grid: L/V dz+l ) TcL/V

(58)

The state variables are just updated with the intermediate approximation:

xz+1 ) xint

(59)

4.2.2. Adaptive Grid Selection. Adaptivity is a general trend, in modern computational science. The objective of adaptive algorithms is to minimize the computational effort to obtain an approximate solution with a prescribed accuracy. For this purpose, adaptive algorithms search for problem formulations (discretizations) with the lowest complexity, for which a prescribed accuracy requirement can just be met. To obtain this optimal discretization, an algorithm that uses a high resolution only in those domains where this high resolution is actually needed has to be developed. In algorithm 1, the discretization has been done on the same grid Λr for each iteration loop. Obviously, uniformity of the grid may detail the grid in domains where the current approximation is already quite good while it might give an unsatisfactory result for other domains with poor approximation. To spare unnecessary effort introduced by redundant refinement, it is possible to change the grid for each step 1 during the iterations. Although an adaptive scheme could certainly be realized with one single adaptively generated set of basis functions Λrz for the discretization of each continuously formulated equation fi, it would leave a large potential for adaptation unused. This potential stems from the fact that it may be highly desirable to choose different sets of basis functions corresponding to different balance volumes in a model. For example, the physical representation can be different in the top and in the bottom of a distillation column, reflecting the different compositions. Generally, each function fi could be associated with its own set of projection functions Ψri . To allow this more flexible choice of approximating projection functions, every appearance of Ψr and Λr in r with algorithm 1 has to be formally replaced by Ψz,i r the corresponding index pair set Λz,i. As will be seen in algorithm 4, the full flexibility will not be exploited in the case of the flash example. However, this flexibility becomes vital if multistage processes such as distillation columns will be considered in future research. Provided this formal replacement has been performed, step 3 in algorithm 1 can be extended to check the error of the current iterate and follow an appropriate adaptation strategy. Each of the steps in the following algorithm has to be performed for each function fi, i ) 1, ..., n. Algorithm 3: Adaptive Grid Selection for the Basic Multigrid Algorithm. Step 3.1. To obtain a measure for the error of the intermediate approximation, a residual evaluation int int rint i (ξ) ) fi(Fz (ξ),xz ,p)

(60)

on the full grid discretization is employed. According to the norm equivalence of the Haar basis, the L2 norm

i )

||rint i (ξ)||L2

(61)

int i ) ||〈fi(Fint z (ξ),xz ,p), Ψ〉||l2 ) ||δdi||l2

(62)

The residual vector δdi contains the residuals for each of the discretized equations:

δdi ) (..., (δdi,(j,k)), ...)T, (j, k) ∈ Λ

(63)

For the flash example, the values of i are measures for the approximation quality of the material balance and the phase equilibrium equations. Step 3.2. Having this measure for the error available, it can be compared against a prescribed tolerance tol by

sii < tol

(64)

The weighting factors si are introduced to level out problem-dependent scaling of the functions fi. If the inequality (64) holds for all i, terminate the calculation with dint and xint as the converged solution; else, z z continue. Step 3.3. The adaptation step consists of several substeps. and Substep 3.3.1. First the new adapted grid Λr,ad i the number of basis functions m in this set need to be initialized:

) { }, m ) 0 Λr,ad i

(65)

Because the new adapted set of the basis function always starts with an empty set, the algorithm allows both refinement and coarsening of the grid. Substep 3.3.2. The adaptation strategy is based on the selection of the residual vector component δdi,(j*,k*), which contributed most to the error norm:

(δdi,(j*,k*))2 ) max{[diag(δdi)]‚δdi}

(66)

According to the discretized system in eq 47, all residual components with indices (j, k) from the current discretization set Λr are zero. Therefore, the selected coefficient can only correspond to a basis function that has not been a member of the current set Λr. If the use of the basis function with the indices (j*, k*) in the next iteration is assumed to have only a minor effect on the vector of residuals δdi for all components with indices (j, k) * (j*, k*), the error ji after the next iteration when (j*, k*) is included in the active set can be estimated:

ji )

x||δd ||

i l2

2

- (δdi,(j*,k*))2

(67)

Because the selection of the “strongest contributer” will be repeated, the corresponding residual vector entry will be set to zero:

(δdi,(j*,k*)) :) 0

(68)

Finally, the discretization set and the number of currently used basis functions are updated:

:) Λr,ad ∪ {( j*, k*)} Λr,ad i i

(69)

m :) m + 1

(70)


Substep 3.3.3. If the scaled estimated error siji > tol and m < mmax, go to substep 3.3.2; else, proceed. Note that the number of added basis functions m is limited to mmax even if the error criterion is not fulfilled. This is necessary in order to prevent the selection of a very large number of basis functions resulting in a poor performance. Step 3.4. Perform problem-specific postprocessing of (see algorithm 4). the adapted sets Λr,ad i Step 3.5. Update the current reduced set of projection functions according to the new adapted grid: r Ψz+1,i :) (..., ψj,k, ...)T, ( j, k) ∈ Λr,ad i

(71)

r :) Λir,ad Λz+1,i

(72)

In algorithm 3, the weighting factors si (see eq 64 and step 3.3.3) have not been defined. These specifications are problem-dependent. For the flash problem, the weighting factors are set to the inverse of the feed flow rate and to unity for the material balance and the vapor-liquid equilibrium, respectively:

s1 :)

1 , s2 :) 1 Pfeed

(73)

This ensures that the weighted errors are at the same order of magnitude for both equations. Though one has the flexibility to use different sets of basis functions for discretization of the different equation, in the case of the simple flash model, it is more reasonable to use only one set of basis functions for discretization because the flash is a one-stage process. Therefore, step 3.4 in algorithm 3 can be further detailed. Algorithm 4: Flash-Specific Extension of the Adaptation Procedure. Step 3.4. Merge the sets of basis functions:

:) Λr,ad ∪ Λr,ad Λr,ad 1 1 2

(74)

Λr,ad :) Λr,ad 2 1

(75)

Furthermore, numerical experiments have revealed that the use of not fully supported sets of basis functions can lead to numerical instabilities of the algorithm. A set of basis functions is called fully supported if the coefficients in the active set form a graded tree in the scalelocation plot (see section 5.2). This occurs if for each wavelet in the active set all wavelets on lower scales with overlapping support are also active. If the sets do not already have this property, it is enforced Λr,ad i by adding the needed basis functions on lower scales. This extension would look different if, for example, a distillation column would be considered. There it would be reasonable to use the same grid for all equations corresponding to one tray. The grids on different trays would, of course, be different. 5. Results Some numerical results will be presented considering the solution of the flash example given in eqs 2-6. Though tailored linear and nonlinear solvers can be expected to enhance the performance of the algorithm, standard software packages have been used in this work. The linear system of eqs 55 and 56 has been

Figure 4. Feed distribution of the two example mixtures.

solved by the package MA48 from the Harwell Subroutine Library (AEA Technologies). A damped NewtonRaphson method implemented in the software package NLEQls30 was used to solve the nonlinear system given in eq 47. The first part of the discussion will show the performance of the algorithm by moving stepwise through the algorithm, while the second part will report results on the effectiveness of the proposed approach. However, before discussion of the results, the following section briefly introduces the used feed mixture and the property calculation methods. 5.1. Mixture Characterization and Property Calculation. To show that the proposed method is not limited by the available experimental characterization or to a continuous thermodynamic framework, two extreme cases of characterization are considered. The minimum characterization necessary to derive compositional description and property calculation methods is given by a distillation curve and the bulk specific gravity. Though this may hardly be satisfactory for a predictive simulation, it represents the minimal experimental information a reasonable simulation can be built on. The most detailed information one can expect from an experimental analysis is the composition of the mixture in terms of all identified real chemical components. For each of the two extreme characterization types, one example mixture is considered. With modifications in the characterization procedure, all other analytical data whose degree of detail may lie between should also be usable in the adaptive multigrid framework. A data set for the minimum characterization was taken from Miquel et al.31 further denoted by MiCa1. They give a TBP curve with 11 data points spanning a boiling point interval from 270 to 900 K. The data set with the detailed composition analysis have been provided by Vacher32 (Institut Franc¸ ais du Pe´trole). This set, in the following denoted by IFP2, contains the detailed GC analysis of a commercial gasoline. To derive a compositional description for the mixture from the minimum characterization, a method proposed by Miquel and Castells33,34 has been extended to allow a continuous composition representation. A detailed description of this extension has been reported by Briesen.22 The two resulting feed distributions are given in Figure 4. Obviously, the TBP characterization results in a smooth distribution while the detailed GC analysis leads to a highly irregular distribution. In both distribution functions, the components are arranged in order of increasing boiling temperature. Because the examples will focus on thermal separation

2344


Figure 5. Intermediate coarse resolution solution Fint 0 (ξ) after the first iteration compared to the high-resolution reference solution.

processes, the vapor-liquid equilibrium has the main influence on the system behavior. The boiling temperature, therefore, qualifies as a characterization variable because it can be seen as an indicator for the phase behavior of a component. However, this particular order of the components is not required. Generally, the components should be ordered according to similarity in the behavior decisive for the particular process. Note that the algorithm works on correction terms rather than on the actual representation. Thus, the ordering of the components in a way that the distribution can be represented on a minimal grid is no issue. The actual representation of the distribution will always be obtained on the full grid. For the property calculations, namely, the phase equilibrium constant, simple correlations have been employed. The phase equilibrium constant is calculated on the basis of the Clausius-Clapeyron equation in combination with Trouton’s rule, resulting in a dependence only on the boiling temperature of the components.35 For a more comprehensive description of the used property calculation methods and the example mixture, we refer to work by Briesen.22 5.2. Algorithm in Operation. At step 0 of algorithm 1, no information on the shape of the unknown vaporand liquid-phase distribution is available. Therefore, a uniform grid with eight basis functions is selected. This roughly can be considered as a simulation with eight pseudocomponents. For this resolution, it can already be expected that the result will not satisfy the accuracy requirements. However, because the number of basis functions, of course, has a direct influence on the performance of the first solution step, the number of basis functions should not be chosen too large. In all of the investigated examples of this work, uniform initial grids of eight basis functions have been sufficient to give a reasonable approximation of the unknown distribution in the first iteration. In steps 1 and 2 of algorithm 1, the correction model is set up and solved. The resulting intermediate distributions are given in Figure 5. For reference a highresolution model solution is also given. For the highresolution model, 128 basis functions have been selected, again corresponding to a simulation considering 128 pseudocomponents. Figure 5 shows that the shape of the distribution can be reflected qualitatively rather than quantitatively because of the coarse initial grid. Consequently, these intermediate approximations do not meet the chosen accuracy constraint of tol ) 10-2. Figure 6 shows the residual error at the intermediate solution distributed over ξ. Because the weighted residual equations are also distributed over ξ, they directly reflect the quality of the approximations. In domains

Figure 6. Weighted residual error distribution for the material balance (left) and the phase equilibria (right).

Figure 7. Scale-location plot of the error distribution.

of ξ where the approximation of the distribution is bad, high residual errors occur. A way of visualizing a multiscale representation is the scale-location plot. For the two error distributions in Figure 6, the scale-location plot is given in Figure 7. Each of the patches corresponds to a coefficient δdj,k of the wavelet expansions. The horizontal axis gives the localization of the corresponding wavelet in the domain ξ ∈ [0, 1]. The scale j of the basis function is represented on the vertical axis. The brightness of the patches corresponds to the absolute value of δdj,k. Because the solution of eq 47 demands that the multiscale coefficients of the basis functions that have been used for discretization must be 0, the patches on the three lowest scales are all black. The bright patches again show the localization of the domains, where the intermediate solution does not adequately represent the reference solution. According to the adaptive grid selection strategy given in algorithm 3, the multiscale coefficients δdj,k, represented in the scale-location plots, are the basis of the grid adaptation. The basis functions that have been selected by the strategy are also represented in Figure 7. All selected basis functions are highlighted with a white frame. The given grid is the final grid after the flash-specific postprocessing in algorithm 4. To obtain the interpolation of the current intermediate solution to the full grid, the linearized eqs 55 and 56 are solved, resulting in the new iterate for the distributions. In Figure 8, this new iterate is given and compared to the reference solutions. For the higher boiling components (ξ > 0.25), the accuracy of the new iterate is already very good even though this is the first iteration loop. However, the high accuracy in this case must also be attributed to the simple process model. These iterates and the adaptively determined discretization grid are now used in the next inner loop (cf. Figure 3). Setting up and solving the correction model again in steps 1 and 2 result in the intermediate solution given in Figure 9. The intermediate solution


Figure 8. Interpolated solution F1(ξ) in the initial outer loop (cf. Figure 3).

Figure 9. Solution F2(ξ) after the first outer loop with an adapted grid.

is nearly indistinguishable from the reference solution. This intermediate solution also satisfies the error constraint, so that it is considered as the “final” solution of the problem. 5.3. Efficiency of the Algorithm. Until now the accuracy of the approximation has been only qualitatively judged. The use of the function norm ||‚||L2, which can be viewed as a generalization of the absolute of a scalar error, provides the means to quantify the error of the obtained approximation. Figure 10 shows this norm for a set of simulations with different model complexities employing two different feed compositions (MiCa1 and IFP2). Besides the adaptive calculations, the figure also includes the results for the direct (fixed-grid) simulation in different resolutions. The direct solutions are obtained on a fixed grid with no further error check or adaptation. The solution for the fixed grid with eight basis functions is the result obtained in the initial step depicted in Figure 5. The model complexity in Figure 10 is measured by the number of basis functions used. Because the results of the direct simulations basically are pseudocomponent calculations, the model complexity can also be thought of as the number of pseudocomponents used for the compositional description of the system. Figure 10 reveals that the error for the direct simulations in both systems decreases to zero for the model complexity of 128 because this is considered as the exact solution of the problem. The rate of decrease, however, is different for the different feed compositions. While for the MiCa1 feed the error decreases quite quickly to reasonable values (Figure 10, top), the simulation with the detailed feed specification IFP2 shows large errors in the compositional description even for large model complexities (Figure 10, bottom). This behavior stems from the fact that this feed distribution has a less smooth shape than the MiCa1 distribution. The adaptive multigrid simulations have been performed with a set of different accuracy specifications. In almost all cases, the error norm for the adaptive

Figure 10. Error of the compositional results for a direct calculation on fixed grids and for adaptive multigrid calculations for different model complexities (specifications: feed flow, 100 kmol/h; feed pressure, 1 bar; feed temperature, 300 K; system pressure, 1 bar; splitting factor R ) 0.5; for adaptive runs, tol ) {2, 1, 0.5, 0.2, 0.1, 0.05, 0.02, 0.01}; for fixed grid, number of pseudocomponents, {2, 4, 8, 16, 32, 64, 128}).

calculations is practically indistinguishable from zero, indicating that a nearly exact solution has been found. Only for the MiCa1 feed, larger errors occur corresponding to the run with the loosest accuracy constraint of tol ) 2. For this run the errors are actually identical to the direct simulation errors. The reason for this is that the constraint was to loose to require an adaptation loop. Consequently, the calculation was terminated after the initial step, which of course results in the same error as that for the direct simulation. Because the system size for the adaptive calculation changes with the iteration, the model complexity needs to be defined differently. For the adaptive calculations, the model complexity is the accumulated number of basis functions which have been active at any iteration stage of the adaptation process. This can be interpreted as an overestimation of the model complexity compared to the direct simulation. The adapted complexities, for which the model has been solved, can be considerably lower. However, in the case of the flash example, the model complexity is the same as the actual model size because all of the simulations needed only one adaptation step. Thus, only two outer iteration steps were needed to obtain the converged result. Considering more complex process models, also the complexity of adaptation behavior should significantly increase. To see how the strategy adapts to different model specifications, Figure 11 shows the results for two different specifications of the splitting factor R. The top and bottom of the figure show the vapor- and liquidphase distribution function, respectively. The center depicts the basis functions that have been used in the course of the discretization. For R ) 0.05 (Figure 11, left), most of the fluid is removed from the bottom of the flash. So, the distribution in the liquid product is

2346


Figure 11. Results and adaptation grid for two different flash specifications (R ) 0.05, left; R ) 0.95, right).

very similar to the MiCa1 feed distribution given in the left-hand side of Figure 4. Instead, the vapor phase contains only the very light end components, giving rise to a strongly increasing distribution for small values of ξ. For the specification of R ) 0.95 (Figure 11, right), the opposite behavior can be observed. Most of the fluid is removed from the top. The vapor-phase distribution is similar to the feed, while the liquid phase contains only the heavy end components. More interesting are the adaptively generated sets of discretization basis functions for both cases. The grids show that the physical behavior is reflected in the discretization grid. In the case of R ) 0.05, the heavy end components do not play any role in the phase behavior, so the grid can be quite coarse for values of ξ > 0.5. For R ) 0.95 the light end still plays an important role in the vapor-liquid equilibrium. However, in contrast to the first simulation, now also the heavier components start to gain influence on the phase equilibrium calculation. This fact is directly reflected in the choice of the basis function, which also includes basis functions for ξ > 0.5 on scales larger than 3. 6. Conclusion Based on the improvements in the analytical techniques for petroleum fraction characterization, the

detail in compositional information will strongly increase in the future. The inclusion of this detailed characterization on a molecular level seems to be the only way to fundamentally enhance the predictive capabilities of models for petroleum processes.20,36 The presented method aims to provide a new simulation technology to enable simulations using this detailed characterization without actually implementing a model based on single-component balances. However, because of the multigrid approach, the final solution can again be obtained on the molecular level. Thus, a result in terms of single chemical species is available. However, the method does not rely on the detailed characterization. Also, a “macroscopic” characterization by means of distillation curves is possible. This underlines the flexibility of the approach. The focus of this work has been more on the algorithmic development of the adaptive techniques than on the solution of a technically relevant problem. Therefore, only a simple model has been used to investigate the primary behavior of the algorithm. The algorithm is capable of providing an errorcontrolled solution with adapted model complexity. In comparison to standard pseudocomponent techniques, the method provides a measure for the error associated with the approximate solution by means of a residual norm of the detailed model formulation. Another possible benefit of the adaptive error control method is its ability to solve the model at various degrees of accuracy. This scalable model accuracy allows one to simulate specific parts of the process in greater detail than others. In a design process, this feature may be interesting in order to investigate a large number of possible configurations on a coarse level of detail and use the most promising subset of configurations without remodeling on a high level of detail for further studies. The presented expositions only consider the processing of hydrocarbon mixtures because of the impressive economical importance of the refining industry. However, the framework is rather general and could possibly also improve the simulation of other complex multicomponent mixture processes. The tailoring necessary to guarantee an efficient solution strategy is very problemdependent and can hardly be foreseen. In a one-stage operation like the one presented here, only the adaptivity with respect to different specifications could have been shown. In multistage operations, one can expect the complexity of the calculation to be further reduced according to an adapted discretization grid on each specific stage. Appendix: Wavelet-Galerkin Projection of the Flash Model To illustrate the general form of the wavelet-Galerkin discretization introduced in section 3.2, this appendix will give the equation resulting from the waveletGalerkin projection of the flash model. The flash model contains two continuously formulated equations which need to be discretized, namely, the material balance (2) and the vapor-liquid equilibrium (3). Because these equations generally arise when modeling thermal separation processes, the equations shown can easily be adapted to more complex models such as distillation columns. Material Balance. First, the derived discretization method should be applied to the material balance equation. Because the projection in the Galerkin scheme


actually is an integration procedure, all of the addends in the material balance can be treated separately. Performing the projection for the material balance equation according to eq 37 for each term leads to the remarkably simple result

0 ) PfeeddF,r - VdV,r - LdL,r

(76)

The simplicity of this vector equation is a result of the linear structure of the material balance with respect to the distribution functions and of the orthonormality of the wavelet basis. Note that in the multigrid framework the actual unknowns in this equation are not the distribution coefficients but the correction terms (see eq 46). Because the single-scale basis is also orthonormal, it leads to a very similar result. Evaluating eq 42 for the material balance results in F

V

0 ) Pfeedc - Vc - Lc

L

(77)

Vapor-Liquid Equilibrium. The orthonormality of the two bases also allows a simple projection of the first term in the vapor-liquid equilibrium equation for both, the single-scale and multiscale, representations:

∫

1 V F (ξ)‚ψj,k 0

V dξ ) dj,k , (j, k) ∈ Λr

∫01FV(ξ)‚φj,l dξ ) cj,lV,

(j, l) ∈ Γ

(78) (79)

For the product of the liquid distribution FL(ξ) and the phase equilibrium distribution K(‚,ξ), only the discretization in the (local) single-scale representation is just as simple. The disjunct support of the single-scale basis allows direct multiplication of the single-scale coefficients according to

∫01{K(‚,ξ) FL(ξ)}‚φj,l dξ 1 ) ∫0 {[(cK)T‚Φ]‚[(cL)T‚Φ]}‚φj,l dξ ) 2 j/2

∫01{([diag cK]‚cL)T‚Φ}‚φj,l dξ

88). This procedure is summarized as follows:

∫01{K(‚,ξ) FL(ξ)}‚ψj,k dξ 1 ) ∫0 {[(cK)T‚Φ]‚[(dL)T‚Ψ]}‚ψj,k dξ

(84)

∫01{[(cK)T‚Φ]‚[(cL)T‚Φ]}‚ψj,k dξ

(85)

FWT

)

) 2 j/2 FWT-1

) 2 j/2

K L ) 2 j/2cj,l cj,l

(82)

vle , (j, l) ∈ Γr ) 2 j/2cj,k

(83)

∫01{(dvle)T‚Ψ}‚ψj,k dξ

vle , ( j, k) ∈ Λr ) 2 j/2dj,k

(86) (87) (88)

Closure Condition. As introduced in eq 8, the function g may also depend on the state functions. However, this dependence is not a direct dependence on the distribution functions but on a scalar value resulting from the operator O applied to the distributions. An example of this is the closure condition in eq 5. If the particular integral operator

O[FV] )

∫01FV(ξ) dξ

(89)

is introduced, the closure condition can easily be represented by the general equation (8). Especially for the multiscale representation of the distribution function, this integration gives a simple result: V ∫01FV(ξ) dξ ) d-1,0

(90)

So, the operator maps the continuous distribution function to one particular distribution coefficient, which gives the closure condition V -1 0 ) d-1,0

(80) (81)

∫01{(cvle)T‚Φ}‚ψj,k dξ

(91)

for the multiscale representation. For the local single-scale representation, the equivalence to the pseudocomponent approach becomes apparent because the closure condition only requires the sum of the weighted coefficients to be equal to 1:

∑ 2-j/2cl,kV ) 1

(92)

(j,l)∈ Γ

The evaluation of the phase equilibrium function K is simpler in the single-scale representation because standard discrete physical property routines may be used in this case. Instead, the multiscale discretization is more elaborate because there is a product of the two distributions K(‚,ξ) and FL(ξ) that needs to be projected. Because the wavelet functions do have overlapping supports, the coefficients cannot be directly multiplied in order to obtain the distribution coefficients of the product. In case the phase equilibrium function K is available in the single-scale representation (eq 84), we can temporarily switch to the single-scale space (eq 85) by the FWT algorithm in order to perform the direct multiplication. After evaluation of the product, the single-scale coefficients of the product (eq 86) can be retransformed into the wavelet space (eq 87), where the projection on the wavelets can be easily performed (eq

Literature Cited (1) Altgelt, K.; Boduszynski, M. Composition and Analysis of Heavy Petroleum Fractions; Dekker: New York, 1994. (2) Krambeck, F. Continuous mixtures in fluid catalytic cracking and extensions. In Chemical Reactions in Complex Mixtures; Sapre, A. V., Krambeck, F., Eds.; Van Nostrand Reinhold: New York, 1991; pp 42-59. (3) Carrier, B.; Rogaiski, M.; Peneloux, A. Choice of pseudocomponents for flash calculations of petroleum fluids. In C7+ Fraction Characterization; Mansoori, G. A., Chorn, L. G., Eds.; Taylor & Francis: New York, 1989; pp 123-136. (4) Neau, E.; Jaubert, J.-N.; Rogalski, M. Characterization of heavy oils. Ind. Eng. Chem. Res. 1993, 32 (No. 6), 1196-1203. (5) Ritzsch, M.; Kehlen, H. Continuous thermodynamics of complex mixtures. Fluid Phase Equilib. 1983, 14, 225-234. (6) von Watzdorf, R. Adaptive Multiskalenverfahren fur die Modelreduktion von thermischen Viel-stofftrennprozessen; VDIFortschritt-Bericht 561; VDI-Verlag: Du¨sseldorf, Germany, 1998; Reihe 3.

2348


(7) Hariu, O.; Sage, R. Crude split figured by computer. Hydrocarbon Process. 1969, Apr. (8) Bowman, J. Distillation of an indefinite number of components. Ind. Eng. Chem. 1949, 41 (No. 9), 2004-2007. (9) Cotterman, R.; Bender, R.; Prausnitz, J. Phase equilibria for mixtures containing very many components. Development and application of continuous thermodynamics for chemical process design. Ind. Eng. Chem., Process Des. Dev. 1985, 24 (No. 1), 194203. (10) Cotterman, R.; Prausnitz, J. Flash calculations for continuous or semicontinuous mixtures using an equation of state. Ind. Eng. Chem., Process Des. Dev. 1985, 24 (No. 2), 434-443. (11) Ratzsch, M.; Kehlen, H.; Schumann, J. Computer simulation of complex multicomponent hydrocarbon distillation by continuous thermodynamics. Fluid Phase Equilib. 1989, 51, 133-146. (12) Daubechies, I. Ten Lectures on Wavelets; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1992. (13) Dahmen, W. Wavelet and multiscale methods for operator equations. Acta Numer. 1997, 6, 55-228. (14) Carrier, J.; Stephanopoulos, G. Wavelet-based modulation in control-relevant process identification. AIChE J. 1998, 44 (No. 2), 341-360. (15) Flehmig, F.; von Watzdorf, R.; Marquardt, W. Identification of trends in process measurements using the wavelet transformation. Comput. Chem. Eng. Suppl. 1998, 24, S491-S496. (16) Sweldens, W.; Schroder, P. Building your own wavelets at home. Wavelets in Computer Graphics; ACM SIGGRAPH Course notes; ACM: 1996; pp 15-87. (17) von Watzdorf, R.; Urban, K.; Dahmen, W.; Marquardt, W. A wavelet-Galerkin method applied to separation processes. In Scientific Computing in Chemical Engineering; Keil, F., Mackens, W., Vofi, H., Werther, J., Eds.; Springer: Berlin, 1996; pp 246252. (18) Daubechies, I. Orthonormal bases of compactly supported wavelets. Comm. Pure Appl. Math. 1988, 41, 906-966. (19) Haar, A. Toward the theory of function systems. (Zur Theorie der orthogonalen Funktionen-Systeme). Math. Ann. 1910, 69, 331-371. (20) Ramage, M. The petroleum industry of the 21st century: The enabling role of technical computing. FOCAPO 1998, 1-18. (21) Hu, S.; Towler, G.; Zhu, X. Combine molecular modeling with optimization to stretch refinery operation. Ind. Eng. Chem. Res. 2002, 41 (No. 4), 825-841. (22) Briesen, H. Adaptive Composition Representation for the Simulation and Optimization of Multicomponent-Mixture Pro-

cesses; VDI-Fortschritt-Bericht 740; VDI-Verlag: Du¨sseldorf, Germany, 2002; Reihe 3. (23) Qian, S.; Weiss, J. Wavelets and the numerical solution of partial differential equations. J. Comput. Phys. 1993, 106, 155175. (24) Dahmen, W.; Kunoth, A.; Urban, K. A wavelet-Galerkin method for the Stokes-equations. Computing 1996, 56 (No. 3), 259-302. (25) Chen, M.-Q.; Hwang, C.; Shih, Y.-P. A wavelet-Galerkin method for solving population balance equations. Comput. Chem. Eng. 1996, 20 (No. 2), 131-145. (26) Liu, Y.; Cameron, I. A new wavelet-based method for the solution of the population balance equation. Chem. Eng. Sci. 2001, 56, 5283-5294. (27) von Watzdorf, R.; Marquardt, W. Fully adaptive model size reduction for multicomponent separation problems. Comput. Chem. Eng. Suppl. 1997, 21, S811-S816. (28) Wesseling, P. An Introduction to Multigrid Methods; Wiley & Sons: Chichester, U.K., 1991. (29) Briggs, W.; Henson, V.; McCormick, S. A Multigrid Tutorial; SIAM: Philadelphia, PA, 2000. (30) Nowak, U.; Weinmann, L. A family of Newton codes for systems of highly nonlinear equations; Technical Report TR-9110; Konrad-Zuse-Zentrum fu¨r Informationstechnik: Berlin, 1991. (31) Miquel, J.; Hernandez, J.; Castells, F. A new method for petroleum fraction, and crude oil characterization. SPE Reservoir Eng. 1992, 265-271. (32) Vacher, P. (Institut Franc¸ ais du Pe´trole). Personal communication, 1998. (33) Miquel, J.; Castells, F. Easy characterization of petroleum fractions (part 1). Hydrocarbon Process. 1993, 101-105. (34) Miquel, J.; Castells, F. Easy characterization of petroleum fractions (part 2). Hydrocarbon Process. 1994, 99-103. (35) Vakili-Nezhaad, G.; Modarress, H.; Mansoori, G. Continuous thermodynamics of petroleum fluids fractions. Chem. Eng. Process. 2001, 40, 431-435. (36) Katzer, J.; Ramage, M.; Sapre, A. Petroleum refining: Poised for profound changes. Chem. Eng. Process. 2000, 41-51.

Received for review August 8, 2002 Revised manuscript received February 12, 2003 Accepted February 28, 2003 IE0206150

An Adaptive Multigrid Method for Steady-State Simulation of

Recommend Documents