Evaluation of Weighted Residual Methods for the ... - ACS Publications

Article pubs.acs.org/IECR

Evaluation of Weighted Residual Methods for the Solution of a Population Balance Model Describing Bubbly Flows: The LeastSquares, Galerkin, Tau, and Orthogonal Collocation Methods Jannike Solsvik* and Hugo A. Jakobsen* Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU), N-7491 Trondheim, Norway ABSTRACT: In dispersed gas−liquid flows, the bubble size distribution plays an important role in the phase structure and interphase forces, which, in turn, determine the multiphase hydrodynamic behaviors, including the spatial profiles of the gas fraction, gas and liquid velocities, and mixing and mass-transfer behaviors. Thus, fluid particle coalescence and breakage phenomena are important for optimal operation of many industrial process units like the bubble column reactors. The population balance equation (PBE) is considered a concept for describing the evolution of populations of countable entities such as the bubbles in the bubble column. In recent studies, the least-squares method has been adopted for the solution of population balance (PB) problems. A favorable property of the weighted residual methods such as the least-squares technique is that the solution of the density function itself can be obtained from the fundamental PBE formulation, that is, not moment formulations. Hence, in this framework, the inner coordinate space is resolved. Although the interest in the least-squares method is a consequence of some favorable properties, the method should be compared to other techniques in the weighted residual family for the solution of integro−differential equations such as the PBE. For this reason, in this study, the performance of the leastsquares method is compared to the orthogonal collocation, Galerkin and tau methods for the solutions of a PB problem describing bubbly flows. On the basis of the present PB model and simulation results, the simple orthogonal collocation method is considered a good candidate. The orthogonal collocation method holds the simplest algebra and, further, is computationally efficient and gives accurate solutions. In particular, the least-squares method suffers from relatively higher numerical error.

■

INTRODUCTION In the simplest configuration the bubble column is a vertical cylinder with the gas entering at the bottom through a gas distributor and with the liquid phase supplied in batch form.1,2 Despite the simplicity of bubble columns and their widespread use for reaction and separation in the process industries, our understanding of the complex flows in these vessels is still very limited. The dispersed phase of the bubble column plays a major role in determining the hydrodynamic behavior of the gas−liquid system. The complexity of the hydrodynamics relates to the evolution of phenomena such as breakage, coalescence, growth, and convective transport of the bubbles. Such events affect the bubble size distribution in the bubble columns; and consequently, set the interfacial momentum, heat and mass transfer fluxes through the contact area. Thus, the bubble size distribution is important for the operation of bubble columns in order to obtain optimal process performance. The complexity of the fluid hydrodynamics and the interactions between transport phenomena and chemical reaction kinetics results in a challenging modeling problem.3,4 Figure 1 sketches a bubble column indicating the evolution of the bubble size distribution. There are different approaches to model multiphase flow systems depending on the physical phenomena of interest and the characteristics of the problem. In the Lagrangian particle tracking approach, the motion of the continuous fluid is described by the Navier−Stokes equations in an Eulerian frame while the particles are tracked through the flow field by solving the equation of motion for each particle, or at least for a number of representative particles. For large systems such as © 2013 American Chemical Society

industrial reactors and when the interaction between the particles becomes important (i.e., dense systems), the Lagrangian approach becomes more complex and computationally expensive. For problems involving complex interactions between the dispersed phase particles, the Eulerian two-fluid models are often used instead of the Lagrangian approach. The use of the two-fluid modeling approach can save a significant amount of computational time compared to the Lagrangian approach.2 In the two-fluid models, the dispersed phase is treated as a continuous medium, analogous to the continuous carrying phase. The two continuous phases assume to coexist in the same physical space and may thus interact. Each phase is governed by its own conservation equations and may interact through phase interaction terms that specify the mass momentum and heat exchanges between the phases. Nevertheless, constitutive equations are required due to the averaging procedure involved in the two-fluid models. Unfortunately, details of the behavior and the characteristics of the interfaces are lost with the two-fluid model formulation. In particular, when the standard two-fluid model is applied to bubbly flows, the dispersed phase is modeled based on an average characteristic bubble size which limits the complexity of the system. A possible alternative is to extend the two-fluid model to a combined multifluid−PBE model. In this approach the dispersed phase is divided into size groups and each group is Received: Revised: Accepted: Published: 15988

June 28, 2013 September 21, 2013 October 16, 2013 October 16, 2013 dx.doi.org/10.1021/ie402033b | Ind. Eng. Chem. Res. 2013, 52, 15988−16013

Industrial & Engineering Chemistry Research

Article

assurance,25 crystallization,26 granulation,27 combustion,28 biological systems,29 etc. Important PBE research activities involve, for example, experimental characterization,30 numerical method studies,31−33 multidimensionality in the internal coordinates,34 and parameter identification.35,36 Some recent PBE research activities of gas−liquid and liquid−liquid flows are summarized in the following. The constitutive equations for breakage and coalescence are crucial in PB modeling. Recent reviews of these models are provided by Solsvik et al.30 and Sajjadi et al.37 Further, simulation analyses of breakage models have recently been performed by Borka and Jakobsen38,39 and Solsvik et al.40 The development of new models for breakage and coalescence is regarded as a significant contribution to the progress of the PB modeling framework.41,42 An important step in the formulation of a PB model for a given process is estimation of values for model parameters. The accuracy of the estimates for these parameters is crucial to the applicability of the formulated model for process design, control, and optimization. Parameter identification studies have been performed by, for example, Raikar et al.,43 Mallikarjunan et al.,44 Jildeh et al.,45 and Singh et al.46 for liquid−liquid systems and by Laakkonen et al.47,48 for gas−liquid systems. Experimental investigation on drop-size distribution for low and high viscous liquid−liquid systems have been performed by Becker et al.49,50 Lucas and co-workers51−57 have provided a number of bubble size distribution data for gas−liquid flows. Unfortunately, not much work have been performed on single drop experiments.58−60 Multivariate PB models include more than one internal coordinate and thus provide the possibility to predict the distributions of different properties of the entities such as size, composition, and temperature.61,62 The combination of fluid dynamics and PBE for the simulation of polydispersed gas−liquid flows has received much attention.62−70 In recent works, the high order least-squares method has been adopted as the solution technique for PB problems.5,71−73 When the PBE is formulated in terms of a density function that depends on the inner coordinate(s), a continuous solution can be achieved by use of a weighted residual method like the leastsquares technique.74 Thus, in this framework, the solution of the density function itself is obtained instead of only a few moments of the density function. As the density function is directly obtained the problem of reconstruction of the density function is avoided.10,14,75 In the work of Nayak et al.74 a Boltzmann-like equation was used to describe the dispersed phase behavior. Distinct from the standard derivation of the moment equations, which are obtained by averaging the Boltzmann-like equation over the velocity space and the internal coordinates in one single operation,76−78 the model formulation presented by Nayak et al.74 split the averaging process into two different operators. The two operators are the averaging in the velocity space and the averaging in the inner coordinate space. The main advantage of the novel derivation by Nayak et al.74 is that one obtains an intermediate result denoting the momentum equations for mass and momentum conservation in terms of the set of internal coordinates. These Maxwellian averaged equations may be employed to resolve the inner coordinate space physics provided that a sufficient numerical method is available. Using the weighted residual approach, like the least-squares method, the inner coordinate space physics may be resolved, and from the distribution function one may compute any of the moments of the inner coordinates in a post-processing procedure. The hypothesis of the suggested model formulation by Nayak et al.74 is that by

Figure 1. Bubble column. Because of different interaction processes, the distribution of the bubble size, ξ, evolves.

considered as a different phase. Thus, more details of the dispersed phase can be achieved by complementing the twofluid model with the PBE.2,5 In the mathematical modeling framework, the PBE is considered a concept for describing the evolution of populations of countable entities such as the bubbles in a bubble column.2,5−11 When the PB modeling technique is applied to bubbly flows, the dispersed gas phase is treated as a population of bubbles distributed not only in physical space but also in an abstract property space. In the PB terminology, the external coordinates refer to the spatial coordinates, and the property coordinates, for example, bubble size, chemical composition, etc., are considered as internal coordinates. The quantity of basic interest in PB modeling is the statistical density distribution function, or density function for short, which represents the behavior of the population of bubbles. The evolution of the density function must take into account the different processes that control the population of the bubbles in the vessel, such as breakage, coalescence, growth, and convective transport. Thus the PBE provides a statistical description of the dispersed phase where the density function may be denoted by f(r,ξ,t) where r is the spatial position vector, ξ is the property of interest of the dispersed phase, and t is the time. Moreover, f(r,ξ,t) dξ represents the average number of particles per unit volume around the point r (r,r + dr) in the time t, with the property between ξ and ξ + dξ. The resulting PBE is characterized mathematically as a nonlinear integropartial differential equation which must be solved by a suitable numerical method.6,12−20 The PBE appears in a variety of applications: aerosol,21 nuclear reactor flow,22 phase separation,23 emulsification,24 flow 15989

dx.doi.org/10.1021/ie402033b | Ind. Eng. Chem. Res. 2013, 52, 15988−16013


Article

Simplification of the Population Balance Equation. In this study, the following simplifications are introduced into the PBE 1: (i) steady state, (ii) cross-sectional averaging giving one independent spatial coordinate z, and (iii) prescribed velocity field of the dispersed phase, that is, vr(r,ξ,t) = vz = constant. In addition, the constitutive eq 9 is used to describe the velocity transporting the particles in the property space, vξ. The LHS of the PBE 1 is thus simplified substituting the constitutive eq 9. With the discussed simplifications, the PBE is presented as

investigating the intermediate equations, the microscopic physics like the particle breakage and coalescence processes may be better understood, and thus it might be possible to derive more reliable and predictive models for these processes. Hence, in this view, the least-squares method may have potentials for detailed studies of the hydrodynamic flows within the bubble columns using a combined multifluid−PBE model. Although much effort has been placed on the least-squares method for the solution of PB problems,5,16,38,39,57,71−74,79−88 other methods in the weighted residual framework are available. The density function itself can also be achieved by use of the orthogonal collocation, Galerkin, and tau methods if the PBE is formulated in terms of a density function that depends on the inner coordinate(s). The orthogonal collocation method has become popular for the solution to chemical reactor engineering problems, such as studies of multicomponent mass diffusion and intraparticle diffusional limitations89−92 and fixed packed bed reactors.93,94 On the other hand, the Galerkin, tau, and least-squares method are relatively seldom observed as the solution technique for chemical reactor problems. However, some applications of the Galerkin method is given by, for example, Nigam,95 Roussos et al.,96 and Shafi et al.97 The performance of the least-squares method for the solution of PB problems should be compared to other methods in the weighted residual framework. For this reason, in this study, methods in the weighted residual family, that is, the orthogonal collocation, Galerkin, tau, and least-squares, are adopted for the solution of a bubble column model. The bubbly flow problem considered is a cross-sectional averaged steadystate model with a prescribed gas-phase velocity. The numerical solution techniques are evaluated on the basis of implementation issues, computational time, residual measures, and condition numbers.

vz

vz ∂ρd (z) = :(ξ , z) 3ρd (z) ∂z

(2)

:(ξ , z) = − b(ξ)fd , m (ξ , z) + V (ξ )

∫ξ

ξmax

h(ξ , ζ )b(ζ )

fd , m (ζ , z)

dζ V (ζ ) 3 fd , m (ζ , z) − ξ 3)1/3 (ξmax dζ − fd , m (ξ , z) c(ξ , ζ ) ξmin ρd (z)V (ζ )

∫

ξ 2V (ξ) + 2 ×

∫ξ

3 1/3 (ξ 3− ξmin )

c([ξ 3 − ζ 3]1/3 , ζ ) [ξ 3 − ζ 3]2/3

min

fd , m (ζ , z) fd , m ([ξ 3 − ζ 3]1/3 , z) ρd (z)V (ζ )

V (ξ ) − V (ζ )

dζ (3)

Constitutive Equations. A set of constitutive equations are required in order to obtain a closed PB problem. The breakage frequency function and the daughter size redistribution function by Coulaloglou and Tavlarides99 together with the coalescence kernel of Prince and Blanch100 are adopted. The breakage frequency function yields99 b(ξ) =

+ ∇r ·[vr(r, ξ , t )fd , m (ξ , r, t )]

⎤ fd , m (ξ , r, t ) ⎡ ∂ρ (r, t ) ⎢ d − + vr(ξ , r, t ) ·∇r ρd (r, t )⎥ ρd (r, t ) ⎣ ∂t ⎦ ∂ + [vξ(r, ξ , t )fd , m (ξ , r, t )] ∂ξ 3 − vξ(r, ξ , t )fd , m (ξ , r, t ) ξ = −:BD(ξ , r , t )

ξvz ∂ρd (z) ∂fd , m (ξ , z) 3ρd (z) ∂z ∂ξ

with the birth and death of particles due to breakage and coalescence given like88

THE POPULATION BALANCE EQUATION On the basis of the work of Nayak et al.,74 a transient PBE holding convective, breakage, and coalescence terms can be written as follows:10,38,39,88,98 ∂t

∂z

−

− fd , m (ξ , z)

■

∂fd , m (ξ , r, t )

∂fd , m (ξ , z)

k1ε1/3 ξ 2/3

⎡ ⎤ σk exp⎢ − 2/32 5/3 ⎥ ⎢⎣ ρl ε ξ ⎥⎦

(4)

where k1 and k2 are empirical parameters that depend on the system properties. The daughter size redistribution function is given as99 h(ξ , ζ ) = 2P(ξ , ζ )

( π2 ξ 2) exp⎛⎜−4.5 [2V (ξ) − V (ζ)]2 ⎞⎟

2.4 =2 (1)

V (ζ )

⎝

V (ζ )2

⎠

(5)

The coalescence closure is given as100

where fd,m(ξ,r,t) is the mass density function which describes the population of bubbles distributed according to their diameter. Two velocity definitions are used in the PBE 1: (i) the velocity denoted by vr(r,ξ,t) represents the velocity field that transports the bubbles in the physical space, and (ii) the velocity denoted by vξ(r,ξ,t) can be interpreted as the velocity field that transports the particles in the property space. For example, in gas−liquid systems, a bubble growth velocity may arise due to the hydrostatic pressure gradient, or by interphase mass transfer. Birth and death of particles due to breakage and coalescence events are collectively denoted by SBD.

π c(ξ , ζ ) = KC [ξ + ζ ]2 [β[εξ]2/3 + β[εζ ]2/3 ]1/2 4 ⎡ ⎡ r 3(ξ , ζ)ρ ⎤1/2 ⎤ 1/3 ⎡ h0 ⎤ l ⎢ ⎢e ⎥ ε ln ⎥ ⎢ ⎥ ⎣ 16σ ⎦ ⎣ hf ⎦ ⎥ × exp⎢ − ⎢ ⎥ [re(ξ , ζ )]2/3 ⎥ ⎢ ⎣ ⎦

(6)

where KC is a calibration factor and re is defined as 15990


Industrial & Engineering Chemistry Research 1⎡1 1⎤ ⎢ + ⎥ 4⎣ξ ζ⎦

Article

−1

re(ξ , ζ ) =

there are algebraic equations. However, the rigorous way to derive the least-squares method is by following the variational principles. In the weighted residual framework, an analogy to the variational principle is adopted for the implementation of the boundary conditions. The different formulations discussed are outlined in appendix based on the generalized problem formulation 13. Details on the least-squares method are given by, for example, Bochev and Gunzburger,104−106 Jiang,103 Bolton and Thatcher,107 Pontaza,108−110 Pontaza and Reddy,111−113 Proot,102 and Proot and Gerritsma.114,115 Outlines of the least-squares theory provided by the chemical engineering community may be found in, for example, Dorao,5 Zhu,71 Sporleder,73 and Solsvik and Jakobsen.98,116 Galerkin, Tau, and Orthogonal Collocation Methods. For particular problems, it is possible to formulate spectral methods starting out from variational principles. The variational methods provide a means to obtain approximate solutions to differential equations based on minimizing the functional of the equation. For example, the Rayleigh−Ritz method117,118 is applied to approximate the solution of variational problems and problems that reduce to variational problems. In particular, the Rayleigh−Ritz method applies calculus of variation for some energy minimization problems due to the fact that the nature organizes itself to a state of minimum energy. The Galerkin method119 can be considered a broad generalization of the Ritz method because it can be used not only as an approximate solution of variational problems, but also for problems that do not reduce to variational problems. The Galerkin method represents a subclass of the method of weighted residuals and is always applicable because it does not depend on the existence of a variational principle. There is always a Galerkin method that corresponds to the variational method, but the reverse is not true. Thus, in one view, the Galerkin method may be interpreted as either variational principle or regular method of weighted residual. To distinguish the two methods, a system of equations is said to be self-adjoint when we can formulate a variational principle for the corresponding problem. Because of the Ritz method is applicable only for equations that are selfadjoint, the Galerkin method is commonly considered a more universal approach for the construction of a spectral approximation. The variational and Galerkin methods give, however, identical results for a system of self-adjoint eqs 2. Let the general framework of the weighted residual methods be presented as

(7)

The volume of bubbles with diameter size ξ is computed from the volumetric relation of an sphere: V(ξ) =

3 πξ 3 4 ⎡ξ⎤ = π⎢ ⎥ 6 3 ⎣2⎦

(8) 73,101

The bubble growth velocity is estimated as vξ(z , ξ) = −

ξvz ∂ρd (z) 3ρd (z) ∂z

(9)

The hydrostatic pressure in the bubble column is computed as p(z) = p0 + ρl g[zmax − z]

(10)

where the effect of the dispersed gas phase on the hydrostatic pressure is neglected. The dispersed phase density can be estimated from the ideal gas law: ρd (z) =

■

ρ (z ) Mdp(z) p(z ) ⇒ d = RT ρ0, d p0

(11)

SPECTRAL SOLUTION TECHNIQUES The goal is to present the PBE 2 on the algebraic system form: (12) Af = F with A and F defined according to the least-squares, Galerkin, tau, and orthogonal collocation techniques. For this task, a generalized boundary value problem is formulated:

3f (z , ξ) = g (z , ξ) )f (z , ξ) = fΓ (z , ξ)

in Ω on Γ

(13a) (13b)

where f(z,ξ) is the unknown function and g(z,ξ) is a given source function. Moreover, 3 • is a linear operator and ) • is the boundary condition operator. The problem is defined in the domain Ω = [zmin,zmax] × [ξmin,ξmax] with the boundary condition applied on Γ, that is, the boundary part of the domain Ω. It should be noticed that, in a similar manner as for one-dimensional problems, a multidimensional problem can be derived and presented on the conventional form 12 using a global index numbering procedure. Basic mathematical prerequisite tools required in the spectral solution framework are outlined in appendix. In particular, local and global indices, trial expansions, differentiation, and quadrature are discussed. Least-Squares Method. There are two conceptually different options available for the formulation of the spectral solution methods: (i) the variational method and (ii) the weighted residual method. Moreover, there are two ways to implement the boundary conditions: (i) the weak formulation and (ii) the strong formulation.102,103 In the weak form, it is required that the residual of both the model equation and the boundary conditions are satisfied at the boundary domain. On the other hand, in the strong form, some linear combinations of the polynomials to satisfy the boundary conditions are performed, or alternatively, the boundary conditions are enforced as additional constraints. In the latter case of the strong form implementation, the algebraic equation system of the model equations is relaxed by the number of boundary condition constraints to give equal number of unknowns as

⟨wI , 9⟩ = 0

∀ I = 0, 1, ..., 7

(14)

where wI represents a weight function and 9 denotes the residual of the generalized problem 13. The choice of weight function distinguishes between the spectral solution schemes. The weight function is identical in the Galerkin and tau methods (i.e., weight function is set equal to the basis function in the polynomial trial function expansion, see eqs 28−31). The essential difference between the Galerkin and the tau methods is the treatment of the boundary conditions (strong form). In the tau method, the boundary conditions are enforced as additional equations and in order to get a system of equations where the number of unknowns is identical to the number of equations, the equation system has to be relaxed by n residual equations and replaced by the n boundary conditions. In the Galerkin method, some linear combinations of the polynomials that fulfill the boundary conditions must be performed. 15991



Article

The simplicity of the orthogonal collocation method120−125 results from employing the Dirac delta function as the weight function. Further details are provided in the appendix. Algebraic System Form. Considering the generalized problem formulation 13, the definitions of the matrix A and vector F in eq 12 for the different solution schemes are presented as

residual1 = ((Af − F)T (Af − F))1/2

(23)

The residual based on the problem operator form is defined as residual 2 = ((Lf − g)T (Lf − g))1/2

(24a)

or alternatively, for the weak formulation of the least-squares method:

Least-squares weak formulation:

residual3 = ((Lf − g)T (Lf − g))1/2

[A]ij = ⟨3Φj , 3Φ⟩ i + ⟨)Φj , )Φ⟩ i

(15)

[F]i = ⟨g , 3Φ⟩ i + ⟨fΓ , )Φ⟩ i

(16)

+ ((Bf − f Γ)T (Bf − f Γ))1/2

(24b)

It should be noticed that the boundary conditions are incorporated in the residual measures 23 and 24. Moreover, for the orthogonal collocation method the residual measures 23 and 24 coincide (see eqs 103 and 104).

Least-squares strong formulation:

■

[A]ij = ⟨3Φj , 3Φ⟩ i

(17)

[F]i = ⟨g , 3Φ⟩ i

(18)

RESULTS AND DISCUSSION PB problems describing bubbly flows are solved using different methods in the weighted residual family, that is, the least-

[A]ij = ⟨3Φj , Φ⟩ i

(19)

Table 1. Physical and Numerical Parameters

[F]i = ⟨g , Φ⟩ i

(20)

Galerkin and tau (strong form):

Pξ Pz ξmin ξmax Zmax vz σ ρc ρ0,d p0 ε k2 β

Orthogonal collocation (strong form):

[A]ij = 3Φj

(21)

[F]i = g

(22)

where Φ denotes the basis polynomials (see eqs 28−31) and ⟨•,•⟩Ω = ∫ Ω•• dΩ. For eqs 17−22, further treatment is need for the boundary conditions. The general algebra and boundary condition treatment of the solution schemes are outlined in appendix along with details of the implementation issues of the PBE. Method Complexity. The orthogonal collocation method is associated with the simplest algebra and implementation issues among the spectral methods considered in the present study. The simplicity results from the use of the Dirac delta function as the weighting function (eq 102). The integral calculation steps (associated with the inner products or the norm definitions) that are required in the least-squares, tau, and Galerkin methods are avoided. Hence, the algebraic system form and the problem operator form coincide (eqs 103−104). On the other hand, the least-squares method is considered to hold the most complex algebraic theory and thus involved implementation issues. Nevertheless, the Galerkin method requires some linear combinations to fulfill the boundary conditions which require some extra work. The benefit with the Galerkin method is the reduced size of the algebraic equation system that needs to be solved.

m m m m/s N/m kg/m3 kg/m3 Pa m2/s3

50; 70 15 0.4 × 10−3 0.0262 3 0.5 0.072 998 1.188 101500 0.392 0.106 2

squares, Galerkin, tau, and orthogonal collocation methods. The physical and numerical parameters applied in the simulations are listed in Table 1. The solution methods are evaluated on both linear and nonlinear cases; that is, considering source terms due to breakage, only a linear PBE is obtained as the nonlinearity arises from the coalescence terms (eq 2). Moreover, different degrees of influence of the breakage and coalescence terms on the performance of the numerical methods are studied. The influence of the breakage and coalescence terms are tuned with k1 and KC in eqs 4 and 6, respectively. If not otherwise specified, the weak formulation of the boundary conditions is applied for the least-squares implementation. The evaluation of the methods are based on residual measures, condition number, computational time, implementation issues, and the mathematical complexity of the methods. However, the simulation results obtained in the present study are based on particular test cases and the observations made should therefore not be generalized. Model Prediction. The present PBE formulation combined with fluid dynamics have been compared to experimental data by Nayak et al.74 In the present study a simplified velocity field is adopted. Hence the phase momentum equations are not considered for a rigorous description of the bubbly flow. Rather, focus is placed on the numerical method performance for spectral solutions of the PBE. Numerical test cases are defined to both showcase and challenge the numerical methods. However, any PB model of interest should be able to predict

■

RESIDUAL MEASURE DEFINITIONS Two residual measures are defined to evaluate the accuracy of the numerical solutions: (i) a residual measure to consider the fulfillment of the model equations, and (ii) a residual measure that considers the algebraic system that is actually solved in order to obtain the solution function. For example, consider the weak formulation of the least-squares method. Equations 72−74 represent the algebraic system form. Moreover, the algebraic system form is given in terms of the problem matrices (L and B) and source vectors (g and fΓ). The residual based on the algebraic system form is defined as 15992



Article

Figure 2. Breakage dominated case. Comparison of the predicted bubble size distribution and Sauter mean diameter with the experimental data T118 by Duan et al.54

Figure 3. Coalescence dominated case. Comparison of the predicted bubble size distribution and Sauter mean diameter with the experimental data M118 by Duan et al.54

study, the focus in not placed on the physics that the models represent but rather on the numerical aspect solving the PBE with weighted residual methods. Only one set of breakage and coalescence kernels are applied in the numerical method analysis: the breakage model proposed by Coulaloglou99 and the coalescence model presented by Prince.100 Different test cases are defined to evaluate the sensitivity to the breakage and coalescence events on the performance of the numerical solution methods. Three test cases, defined by the breakage frequency parameter k1 (eq 4) and the coalescence calibration factor KC (eq 6), are given for the nonlinear breakage− coalescence PBE: (i) a breakage dominated case, (ii) a coalescence dominated case, and (iii) an intermediate case. Further, two test cases are defined for the linear breakage PBE (i.e., KC = 0): (i) slightly influenced by breakage, and largely influence by breakage. The numerical test cases are summarized in Table 2. The Nonlinear Breakage−Coalescence Population Balance Model. The nonlinearity behavior of the PBE arises due to the coalescence terms (see eq 2). An under-relaxation factor of 0.9 was used in the simulations. Moreover, the nonlinear terms were linearized using the Picard iteration approach.131 The numerical test cases A−C in Table 2 are considered in the sequel. Residual. The residual (algebraic system form, eq 23) versus the number of iterations is shown in Figures 4d, 5d, and 6d for cases B, A, and C in Table 2, respectively. Both the tau and Galerkin methods reach the machine precision. Moreover, there is insignificant differences in the solution obtained with the tau and Galerkin method. On the other hand, the computed residual of the orthogonal collocation and least-squares methods are significantly higher than that obtained with the tau and Galerkin methods. However, the number of iterations

Table 2. Numerical Test Cases case A B C D E

k1 (eq 4) 0, 0, 0, 0, 0,

0202 0067 0101 0101 0403

KC (eq 6)

description

figure

10−2 3 × 10−2 10−2 0 0

nonlinear; breakage dominated nonlinear; coalescence dominated nonlinear; intermediate case linear linear; strongly breakage inuenced

5 4 6 13 14

trends in experimental data. Unfortunately, only a few crosssectional averaged experimental data sets exist for bubbly flow.54,126 For example, Duan et el.54 provide cross-sectional averaged bubble size distributions under complex gas−liquid bubbly flow conditions. Both breakage- and coalescencedominated flow conditions were studied. Figures 2 and 3 hold the bubble-size distribution and Sauter mean diameter for breakage- and coalescence-dominated cases, respectively. A Gaussian distribution was adopted in both cases to specify the bubble size distribution at the bed entrance used. However, numerical stability limits the uniformity (i.e., small value of the variance in the Gauss distribution) of the specified bubble size distribution. Moreover, care must be taken at the lower and upper bubble size limits in order to conserve mass in the system.127 The breakage and coalescence models are central in the PBE. The prediction of the shape of the bubble size distribution is importantly related to the choice of constitutive equation for breakage and coalescence. In general, Figures 2 and 3 reveal that the PB model is able to predict the general trends of the experimental data fairly. Breakage and Coalescence Events. There is a variety of breakage and coalescence models suggested in the literature, for example, the reviews by Liao and Lucas,128,129 Lasheras et al.,130 Sajjadi et al.,37 and Solsvik et al.30 However, in the present 15993



Article

Figure 4. Solution of the PBE holding the convection, growth, breakage and coalescence terms. Case B in Table 2 (coalescence dominated).

to reach convergency differs among cases A−C in Table 2. The highest number of iterations are required for the coalescence dominated case, as is reasonable. The least-squares method is most sensitive to the degree of influence of the breakage and coalescence events. Among cases B, A, and C in Table 2, the least-squares method obtains the largest residual for the coalescence dominated case. Considering the Figures 4d, 5d, and 6d, the Galerkin and tau methods hold the fastest convergency rate, followed by the least-squares method. On the basis of the residual measure definition, eq 24, the Galerkin and

tau method have clearly favorable convergency properties compared to the least-squares and orthogonal collocation methods. With the use of the residual definition, eq 24, the residual is directly related to the fulfillment of the PBE 2 and the boundary conditions. By adopting the residual measure based on the problem operator form, eq 24, completely different residual plots are obtained for cases B, A, and C in Table 2 (Figures 4e, 5e, and 6e) compared to the residual measure based on the algebraic system form 23 (Figures 4d, 5d, and 6d). 15994



Article

Figure 5. Solution of the PBE holding the convection, growth, breakage and coalescence terms. Case A in Table 2 (breakage dominated).

On the basis of the residual definition, eq 24, there is insignificant differences between the orthogonal collocation, tau, and Galerkin methods. It should be recalled that the same results are obtained for the orthogonal collocation method for residual definitions, eq 24, and 23 (see eqs 103 and 104). However, a remarkable higher residual is obtained with the least-squares method compared to the other methods. The convergency properties of the least-squares method may improve if the Picard iteration technique is replaced with, for example, the Newton iteration method.85

Figure 7 presents residual plots obtained with the leastsquares method. Both the strong and weak formulation of the boundary conditions are evaluated and the effect of the growth term is revealed. In particular, the residual plots in Figure 7a are based on the algebraic system form (eq 23) and the residual plots in Figure 7b are based on the problem operator form (eq 24). The growth term has a significant influence on the residual, for both residual definitions 23 and 24. If the growth term is excluded from the PB model, a lower residual can be obtained with the least-squares method. On the other hand, the 15995



Article

Figure 6. Solution of the PBE holding the convection, growth, breakage and coalescence terms. Case C in Table 2.

corresponding residual plots of the tau method in Figure 8 are not influenced by the growth term. Figure 12 shows that the evolution of the bubble size distribution due to convection in the property coordinate can be considered significant. However, the differences in the residual obtained by the strong and weak formulation of the boundary condition in the leastsquares method do not differ for the residual definition 23 (algebraic system form). On the other side, for the residual measure based on the problem operator form, eq 24, the strong formulation of the boundary conditions is favorable if the

growth term is not included. The differences in the solution adopting the weak and strong forms are revealed in Figures 9 and 11. Deviations are observed between the tau method and the weak formulation of the least-squares method at the boundary of the ξ-domain. Except for this observation, no observable deviations are present for the solution of the present PB problem (Figures 9, 10, and 11). Dorao5 presents the convergency plot of a variety of PB problems, in which most of the test cases have an analytical solution. In many cases, high accuracy, that is, close to the 15996



Article

Figure 9. Effect of the growth term. The the tau and LSQ (weak form) methods. Case A in Table 2.

Figure 10. Effect of the growth term. The tau and orthogonal collocation methods. Case A in Table 2.

Figure 7. Influence of the growth term on the residual. Least-squares method with strong and weak formulations of the boundaries. The PBE holds the convection, breakage, and coalescence terms. Case C in Table 2.

study. However, the difference in the residual measure definition by Patruno et al.132 and eq 24a should be noted: residual4 = (⟨Lf − g, Lf − g⟩)1/2

machine precision, can be obtained with the least-squares method. High accuracy is also observed in the convergency plots presented by Zhu et al.85 for a breakage−coalescence PB problem with an analytical solution available. For a more comprehensive and physical realistic problem, that is, a twofluid three-phase model describing the mass and momentum conservation for each phase with the PB framework replacing the mass conservation equation for the dispersed phase, Patruno et al.132 presents a convergency plot obtained with the least-squares method that is relatively far from the machine accuracy and thus support the observations made in the present

(25)

In general, the latter residual measure definition gives a lower residual estimate than eq 24a. For both a linear and a comprehensive nonlinear pellet models,133 discuss the limitation of the least-squares method for strongly diffusion limited processes. Significantly higher residual estimates were obtained for the least-squares method compared to the orthogonal collocation, tau, and Galerkin methods using residual measure definition, eq 24. Hence, similar trends as observed in the present study for the PBE are also observed for the pellet equations. The convergency plots presented by Rout134

Figure 8. Influence of the growth term on the residual. Tau method. The PBE holds the convection, breakage, and coalescence terms. Case C in Table 2. 15997



Article

The least-squares method suffer from significantly higher condition numbers than the other methods. The problem with high condition number for the least-squares method is already well-known, for example, Jiang103 and Solsvik.116 The favorable lower condition numbers obtained using the Galerkin method may be a consequence of the reduced equation system resulting from the linear manipulations related to the boundary conditions (see the section Galerkin, Tau, and Orthogonal Collocation Methods). The Linear Breakage Population Balance Model. The effect of the coalescence mechanism is neglected in the linear PB problem. The linear PB model holds the convection terms in the physical and property spaces, and the breakage birth and death terms. In the linear PB problem, no linearization, and thus no iterations, are required and the solution is obtained solving the algebraic equation system only once. The numerical test cases D and E in Table 2 are considered in the sequel. Many of the same observations in the nonlinear breakage− coalescence PB problem are also observed in the linear PB problem. The simulation results of the linear PB problem are given in Figures 13 and 14. As in the nonlinear breakage−coalescence PB problem, the residual definition adopted is of significant importance for the evaluation of the weighted residual methods. If the residual definition based on the algebraic system form (eq 23) is adopted, the Galerkin and tau methods are favorable (Figures 13d and 14d). On the other hand, if the residual definition based on the problem operator form is employed (eq 24), there is negligible differences between the estimated accuracy of the orthogonal collocation, tau and Galerkin methods, whereas the least-squares method suffer from large residuals (Figures 13e and 14e). In particular, whereas there are minor differences in the estimated residuals of the orthogonal collocation, tau and Galerkin methods for cases D and E in Table 2, the leastsquares method is sensitive to the numerical challenging test case E. The trends in the estimated condition numbers (Figures 13c and 14c) follow that of the nonlinear breakage−coalescence PB problem. Further Work. Zhu et al.85 investigate the potential of the least-squares method with a direct minimization algorithm. The solution scheme is analyzed on a simplified PB problem with an analytical solution available. Due to the promising results by Zhu et al.,85 a comparison of the present simulation results of the conventional least-squares method with the least-squares technique with a direct minimization algorithm would be interesting.

Figure 11. Effect of the growth term. The tau and LSQ (strong form) methods. Case A in Table 2.

Figure 12. Solution of the PBE holding the convection and growth terms.

resulted from the solutions of a comprehensive pellet model are questionable as the identical convergency plot is presented for different chemical processes, that is, the methanol synthesis, sorption-enhanced steam methane reforming process, and the steam methane reforming process, which are quite different in nature, and moreover, both spectral and spectral-element formulations of the least-squares method were used. Computational Cost. The computational cost is an important aspect of a numerical method. The evaluation of the computational costs are presented in terms of the relative computational time. Three different approaches are employed: (i) time required to perform a fixed number of iterations (40 iterations), (ii) time used to reach a residual less than a fixed value (10−9) based on the algebraic system form 23, and (iii) time to reach a residual less than a given value (10−9) based on the problem operator form, eq 24. For a given number of iterations (Figures 6f, 4f, 5f) the orthogonal collocation method is the most computational efficient method, whereas the leastsquares is the most time-consuming method. Favorable fewer computation steps are necessary in the orthogonal collocation method due to the property of the Dirac delta function applied in the inner product. On the other hand, if the computational time is evaluated based on case ii, Figures 6g, 4g, and 5g reveal that the Galerkin and tau methods are favorable because of the faster convergency rate. On the other hand, if the computational time is based on case iii, the orthogonal collocation method is slightly faster than the tau method (Figures 6h, 4h, and 5h). Hence, the conclusion of the most efficient solution method is largely dependent on the convergency criteria. However, it should be noticed that the computational time may vary between different computers so the computed relative computational time should be considered a rough estimate. Condition Number. The condition numbers obtained in cases A−C in Table 2 are presented in Figures 6c, 4c, and 5c.

■

CONCLUSION The prediction of the dispersed phase plays a major role in multiphase chemical reactor engineering. The PBE is an integro−differential equation for describing the evolution of the dispersed phase. Many solution methods are suggested in the literature for PB problems. In particular, the least-squares method has recently been applied to solve PB problems. The interest in this method is due to some particular favorable numerical properties. Nevertheless, with the PBE formulated in its fundamental form, that is, in terms of a density function, the distribution function itself is obtained instead of, for example, a few moments. Besides the least-squares approach, the Galerkin, tau, and orthogonal collocation method are high order techniques that can provide detailed solutions to the fundamental PBE formulation. The least-squares method has 15998



Article

Figure 13. Solution of the PBE holding the convection, growth, and breakage terms. Case D in Table 2.

not yet been compared to other weighted residual methods for the solution of the PBE. The motivation of the present study was therefore to reveal the performance of the least-squares method, relative to the Galerkin, tau, and orthogonal collocation method, for the solution of a PBE. A simplified model for bubbly flow was selected for this investigation. Both the computational time and the residual obtained by the weighted residual method are strongly dependent on the criteria used for the evaluation. In the present study, two residual measure definitions have been evaluated. Moreover, the time used to finish the iteration loop has been evaluated based on the two residual measure definitions but also based on a fixed number of iterations. Different conclusions can be made on the relative performance of the weighted residual methods for the PB problem. In general, evaluating a residual related to the solution of the governing equation(s) is considered advantageous above the algebraic system form. Thus, based on the residual measure definition 4 the orthogonal collocation method obtains negligible differences in the accuracy to that of the tau and Galerkin methods. Moreover, the orthogonal collocation method is the most computationally efficient method. The last-squares method, on the other hand, suffers from significantly higher residual (eq 24) and is the most

computational demanding method considered in the present study. The least-squares method needs to be further investigated to elucidate the remarkable residual behavior. The observations made in the present study may not be generalized to other chemical reactor problems. The leastsquares method is still at early stage in the present discipline and further investigations are needed. Thus, further research is required to elucidate the potential of the least-squares method to chemical reactor problems. In particular, the problem associated with the relatively high residual of the problem operator form should be elucidated.

■

MATHEMATICAL PREREQUISITE TOOLS

Lagrange Polynomials

The interpolating polynomial in the Lagrange form is a linear combination of Lagrange basis polynomials: l Pz(z) =

Pz

∑ f j l jP (z) z

jz = 0

15999

z

z

(26)



Article

Figure 14. Solution of the PBE holding the convection, growth, and breakage terms. Case E in Table 2.

in which lPz(z) is the interpolating polynomial of degree Pz and lPjzz(z) is the Lagrange basis polynomials of degree Pz. The Lagrange basis polynomials may be written in product form as l jPz(z) = z

Pz

∏ iz = 0, iz ≠ jz

a polynomial trial function expansion. For a one-dimensional case, the approximate solution can be given as f (z) ≈ f Pz (z) =

z − z iz z j − z iz z

Pz

∑ α j Φ Pj (z) z

z

z

jz = 0

(28)

where ΦPjzz(z) are the basis or trial functions and αjz are the basis coefficients. When the polynomial trial function expansion is substituted into the problem equation and upon the minimization of the residual function, the unknown coefficients αjz are determined. Normally, multidimensional basis functions can be defined as a product of one-dimensional basis functions. For example, the approximate solution of a two-dimensional problem can be expressed as

(27)

Collocation Points

In connection with the solution approximations, the term collocation points is normally associated with the act of assigning the error or residual at given points. If nodal basis functions are chosen for the implementation of a spectral method, the solution function is approximated at the collocation points. The location of the collocation points in the Ω-domain is generally selected as the roots of one of the orthogonal polynomials in the Jacobi family.135 In particular, the Legendre polynomials, which are part of the class of orthogonal polynomials in the Jacobi family, are employed in the present study.

f (z , ξ ) ≈ f 7 (z , ξ ) =

Pz

Pξ

∑ ∑ α j j Φ Pj (z)Φ Pj (ξ) zξ

jz = 0 jξ = 0

z

ξ

z

ξ

(29)

In general, there are two types of basis functions called the modal and nodal bases. Nodal basis functions are commonly used for the implementation of spectral methods because of the resulting simplicity of the method. The nodal basis functions

Polynomial Trial Function Expansion

Briefly, spectral methods are based on the use of a representation of the solution f throughout the domain Ω via 16000



Article

Similarly, for the two-dimensional solution approximation 31 the differentiations are presented like ∂f 7 (z , ξ) = ∂z

∂f 7 (z , ξ) = ∂ξ

Pz

Pξ

∑∑ jz = 0 jξ = 0

Pz

fj j

zξ

dl jPz(z) z

dz

jz = 0 jξ = 0

(33) P

Pξ

∑∑

P

l jξξ(ξ)

f j j l jPz(z) zξ z

dl jξξ(ξ) dξ

(34)

Integral Approximation

A quadrature rule is an approximation of the integral of a function, usually stated as a weighted sum of function values at specified points within the domain of integration. The Gauss quadrature rules136,137 are applied in the present study. Global Index Procedure

The system of algebraic equations of a multi-dimensional problem can be assembled on the conventional form, eq 12, which simplifies the numerical manipulation. In this approach, the numerical solution approximation is expressed by a global index numbering procedure. Given a two-dimensional problem with coordinates ξ and z, the global index J is defined in term of the local indices jz and jξ. By following the global grid labeling as sketched in Figure 15, the global index is given as

Figure 15. Global grid labeling procedure considering a twodimensional case. A GL−GLL-grid where xI is the global numbering, and ξi and zi denote the two local one-dimensional numberings.

are generally defined in terms of a set of Lagrange basis polynomials, eq 27, through a particular set of node values. With this particular type of basis function, that is, ΦPjzz(z) = lPjzz(z) and ΦPjξξ(z) = lPjξξ(z), the polynomial trial function expansions 28 and 29 are given as, respectively

J = jξ + jz (Pξ + 1)

where J ∈ [0, 7 = (Pξ + 1)(Pz + 1) − 1], 0≤jξ≤Pξ and 0 ≤ jz ≤ Pz. Using the global index procedure, the solution approximation 31 can be rewritten in the following way:

Pz

f (z) ≈ f Pz (z) =

∑ f j l jP (z) z

jz = 0

z

z

f ( z , ξ ) ≈ f 7 (z , ξ )

(30)

Pz

and f ( z , ξ ) ≈ f 7 (z , ξ ) =

Pz

=

Pξ

z

zξ

z

ξ

Pz

(31)

=

Because of the Lagrangian polynomial takes the value one at node zjz zeros at all other nodes the coefficients αjz in eq 28 are identified as exactly the solution function values f jz at the node point zjz (eq 30). Hence the problem of finding f 7 (z , ξ) is transformed into the problem of finding unknown basis coefficients. Similarly, for the two-dimensional case the coefficients αjzjξ take the function values f jzjξ at the node points zjz and ξjξ (eq 31).

jz = 0

ξ

ξ

∑ ∑ f j7j l j7j (z , ξ) zξ

zξ

∑ f J7 lJ7(z , ξ) J=0

(36)

in which P

l jPz(z)l jξξ(ξ) = l j7j (z , ξ) = lJ7(z , ξ) z

zξ

(37)

and

f j7j = f J7

(38)

zξ

To indicate that there is a correlation between the local and global indices the following type of notation is adopted in the Pz sequent: f j7(J )j (J ), lPξ jξ(J)(ξ) and ljz(J)(z).

dl jPz(z)

z

z

dz

z

P

=

For ordinary or partial differential problems, the discrete approximation to the derivatives is obtained by differentiation of the approximated solution presented by a polynomial trial function expansion. For a nodal basis in terms of the Lagrange basis polynomials, the differentiation of the one-dimensional solution approximation 30 is given as

z

zξ

Pξ

jz = 0 jξ = 0

Differentiation

∑ fj

z

ξ

lPjzz(z)

Pz

Pξ

∑ ∑ f j7j l jP (z)l jP (ξ) jz = 0 jξ = 0

∑ ∑ f j j l jP (z)l jP (ξ) jz = 0 jξ = 0

df Pz (z) = dz

(35)

ξ

In the global numbering procedure sketched in Figure 15 the ξ points in the grid can be given in vector form as

(32) 16001



⎡ ξ0 ⎤ ⎢ ⎥ ⎢ ξ1 ⎥ ξ=⎢ ⎥ ⎢⋮ ⎥ ⎢ξ ⎥ ⎣ Pξ ⎦

local

⎡ ξ0 ⎤ ⎡ ξ0 ⎤ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ξ ⎢ ξPξ ⎥ ⎢ Pξ ⎥ ⎥ ⎢ ⎥ ⎢ξ ⎢ ξ0 ⎥ ⎢ Pξ + 1 ⎥ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎢ ⎥ ⎢ ⎥ = and ξ = ⎢ ξPξ ⎥ ⎢ ξ2P + 1⎥ ⎢ ⎥ ⎢ ξ ⎥ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ ⎢ ξ0 ⎥ ⎢ ξ7 − Pξ ⎥ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⎣ ξPξ ⎥⎦ ⎢ ξ ⎣ 7 ⎦global

Article

WIQ = W iξqξ(IQ )Wizqz(IQ )

where Wiξqξ(IQ) and Wizqz(IQ) are the local one-dimensional quadrature weights for the coordinates ξ and z, respectively. Further, the local quadrature points are indicated by the indices iqξ(IQ) and iqz(IQ), whereas the global quadrature points are denoted by IQ. Spectral Method Algebra

The generalized problem formulation eq 13 is applied to outline the algebra of the least-squares, Galerkin, tau and orthogonal collocation methods. For the least-squares method, the algebra is outlined based on both the variational principles and the weighted residual framework. Moreover, both the weak and strong formulation of the boundary conditions in the leastsquares technique are illustrated. The algebra of the Galerkin, tau, and orthogonal collocation methods are outlined based on the weighted residual framework. The Least-Squares Method. The Least-Squares Method Based on Variational Principles (Weak Formulation). The basic idea of the least-squares method is to minimize the integral of the square of the residual over the computational domain. In order to find the approximation in the sense of the least-squares error, the following norm-equivalent least-squares function is defined for the generalized problem formulation 13:

(39)

whereas the z-points in the grid is given as such:

⎡ z0 ⎤ ⎢z ⎥ 1 z = ⎢⎢ ⎥⎥ ⋮ ⎢ ⎥ ⎣ z Pz ⎦

local

⎡z ⎤ ⎡ z0 ⎤ ⎢ 0 ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⋮ ⎢ ⎥ ⎢z ⎥ P ⎢ z0 ⎥ ⎢ ξ ⎥ ⎢ ⎥ ⎢ z P +1 ⎥ z ξ ⎢ 1⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ ⎥ z z ⎢ 1 ⎥ ⎢ 2Pξ + 1 ⎥ ⎢ ⎥ and z = ⎢ z 2 ⎥ = ⎢ z 2(Pξ + 1)⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢z ⎥ ⎢ ⎥ 2 ⎢ ⎥ ⎢ z 3Pξ + 1 ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ z ⎥ ⎢⋮ ⎥ ⎢ Pz ⎥ ⎢ z 7 − P ⎥ ξ ⎢⋮ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⎣ z Pz ⎥⎦ ⎢⎣ z 7 ⎥⎦

1(f ; g , fΓ ) =

f7

1 ∥3f (z , ξ) − g (z , ξ)∥Y2 (Ω) 2 1 + ∥)f (z , ξ) − fΓ (z , ξ)∥Y2 ̂(Γ) 2

(43)

with the norms defined as

global

∥•∥Y2 (Ω) = ⟨•, •⟩Y (Ω) =

∫Ω • • dΩ

(44)

∥•∥Y2 ̂(Γ) = ⟨•, •⟩Y ̂(Γ) =

∫Γ • • dΓ

(45)

(40)

Hence, it follows that the solution vector is expressed as ⎡f ⎤ ⎡f ⎤ ⎢ jz=0 , jξ=0 ⎥ ⎢ J=0 ⎥ ⎢ ⎥ ⎢f ⎥ ⎢ f jz=0 , jξ=1 ⎥ ⎢ J=1 ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢f ⎥ ⎢ fJ=P ⎥ ⎢ jz=0 , jξ=Pξ ⎥ ξ ⎢ ⎥ ⎢ ⎥ ⎢ ⎢ fj ,j ⎥ fJ=P +1 ⎥ ξ ⎢ ⎥ ⎢ z=1 ξ=0 ⎥ ⎢ ⎥ = = ⎢⋮ ⎥ ⋮ ⎢ ⎥ ⎢ ⎥ ⎢ f J = 2P ⎥ ⎢ fj ,j ⎥ ξ ⎢ ⎥ ⎢ z=1 ξ=Pξ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ f J = 7 − Pξ ⎥ ⎢ fj ,j ⎥ ⎢ ⎥ ⎢ z=Pz ξ=0 ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⎣ f J = 7 ⎥⎦ ⎢ fj ,j ⎥ global ⎣ z=Pz ξ=Pξ ⎦local

(42)

where the spaces Y(Ω) and Ŷ(Γ) provide the norms for measuring the residuals. On the basis of variational analysis, the minimization statement is equivalent to the following: Find f(z,ξ) ∈ X(Ω) such that lim

d1(f (z , ξ) + εv(z , ξ); g (z , ξ), fΓ (z , ξ)) dε

ε→ 0

∀ v(z , ξ) ∈ X(Ω)

=0 (46)

where the X(Ω) is the solution space. The space X(Ω) serves as a trial space for the candidate minimizers of the functional. Thus, the solution procedure consists in seeking the solution f in the function space X(Ω). The space X(Ω) is determined by the problem and the mathematical formulation used. The interpretation of the minimization statement 48 is the search for a minimum limit of the least-squares functional 1(f + εv ; g , fΓ ) that is expressed as a function of a perturbed solution function f(z,ξ) defined by a generalized perturbation function v(z,ξ) multiplied by a small amplitude variable ε. Thus, in terms of the functional 43 the variational minimization statement 46 is presented by

(41)

Moreover, for integral approximations the two-dimensional quadrature weights are defined as 16002


Industrial & Engineering Chemistry Research lim

ε→ 0

Article

d ⎡1 ⎢ ∥3[f (z , ξ) + εv(z , ξ)] − g (z , ξ) dε ⎣ 2 ⎤ 1 ∥Y2 (Ω) + ∥)[f (z , ξ) + εv(z , ξ)] − fΓ (z , ξ)∥Y2 ̂(Γ) ⎥ ⎦ 2

X7(Ω) = span{l07(z , ξ), ..., l 77(z , ξ)}

Thus, in this particular subspace, any function f 7 (z , ξ) ∈ X7(Ω) can be expressed by eq 36. Inserting the polynomial trial function expansion in terms of Lagrangian polynomials, lj, that is, approximation 36, into eq 51, and choosing systematically

(47)

=0

Provided that 3 • and ) • are linear operators, the variational minimization statement 47 can be rewritten as lim

ε→ 0

v(z , ξ) = l07(z , ξ), ..., l 77(z , ξ)

d ⎡1 ⎢ ∥3f (z , ξ) + ε 3v(z , ξ) − g (z , ξ) dε ⎣ 2 ⎤ 1 ∥Y2 (Ω) + ∥)f (z , ξ) + ε )v(z , ξ) − fΓ (z , ξ)∥Y2 ̂(Γ) ⎥ ⎦ 2

a solvable system of algebraic equations is achieved. The first algebraic equation is achieved by choosing the first element v(z , ξ) = l07(z , ξ) in the subspace X7(Ω) = span{l07(z , ξ), ..., l 77(z , ξ)}

(48)

=0

as such:

Using the definitions of the norms 44 and 45 and further noting that the order of the differentiation and integration operations may be interchanged as the integral limit Ω is not a function of ε, the latter expression can be given as

7

(( ∑ f J lJ7(z , ξ), v(z , ξ) = l07(z , ξ)) J=0

= -(g (z , ξ), fΓ (z , ξ), v(z , ξ) = l07(z , ξ))

1 d lim (3f (z , ξ) + ε 3v(z , ξ) − g (z , ξ))2 dΩ 2 ε → 0 dε Ω 1 d + lim ()f (z , ξ) + ε )v(z , ξ) − fΓ (z , ξ))2 dΓ 2 ε → 0 dε Γ

∫

which is equivalent to 7

∫

= +

⟨3[ ∑ f J lJ7(z , ξ)], 3[l07(z , ξ)]⟩Y (Ω) J=0

∫Ω ([3f (z , ξ) − g(z , ξ)]3v(z , ξ)) dΩ ∫Γ ([)f (z , ξ) − fΓ (z , ξ)])v(z , ξ)) dΓ = 0

7

+ ⟨)[ ∑ f J lJ7(z , ξ)], )[l07(z , ξ)]⟩Y ̂(Γ) J=0

(49)

= ⟨g (z , ξ), 3[l07(z , ξ)]⟩Y (Ω)

In the bracket form of the norms 44 and 45, the minimization statement 49 is given as

+ ⟨fΓ (z , ξ), )[l07(z , ξ)]⟩Y ̂(Γ)

⟨3f (z , ξ), 3v(z , ξ)⟩Y (Ω) + ⟨)f (z , ξ), Bv(z , ξ)⟩Y ̂(Γ)

(55)

We continue to chose elements v(z , ξ) = lI7 in the space v(z , ξ) = {l07(z , ξ), ..., l 77(z , ξ)}

= ⟨g (z , ξ), 3v(z , ξ)⟩Y (Ω) + ⟨fΓ (z , ξ), Bv(z , ξ)⟩Y ̂(Γ) (50)

and moreover, we utilized the linearity of the operators 3 • and ) •. Thus, the latter equation can be rewritten and presented in the general form:

The latter statement is normally rearranged in the matrix form: ((f (z , ξ), v(z , ξ)) = -(g (z , ξ), fΓ (z , ξ), v(z , ξ))

(54)

(51)

7

∑ f J ⟨3[lJ7(z , ξ)], 3[lI7(z , ξ)]⟩Y (Ω)

with

J=0

((f (z , ξ), v(z , ξ))

7

= ⟨3f (z , ξ), 3v(z , ξ)⟩Y (Ω) + ⟨)f (z , ξ), Bv(z , ξ)⟩Y ̂(Γ)

+

∑ f J ⟨)[lJ7(z , ξ)], )[lI7(z , ξ)]⟩Ŷ(Γ) J=0

(52)

= ⟨g (z , ξ), 3[lI7(z , ξ)]⟩Y (Ω)

-(g (z , ξ), fΓ (z , ξ), v(z , ξ))

+ ⟨fΓ (z , ξ), )[lI7(z , ξ)]⟩Y ̂(Γ)

= ⟨g (z , ξ), 3v(z , ξ)⟩Y (Ω) + ⟨fΓ (z , ξ), Bv(z , ξ)⟩Y ̂(Γ)

(56)

Hence, more compactly, the algebraic equation system is allocated on the form 12 with the matrix A and the vectors f and F defined as

(53)

For spectral methods, the solution f(z,ξ) throughout the domain Ω is represented via a polynomial trial function expansion. The discretization procedure thus consists in finding the approximate solution in a reduced subspace, that is,

[A]IJ = ⟨3[lJ7(z , ξ)], 3[lI7(z , ξ)]⟩Y (Ω) + ⟨)[lJ7(z , ξ)], )[lI7(z , ξ)]⟩Y ̂(Γ)

f 7 (z , ξ) ∈ X7(Ω) ⊂ X(Ω)

(57)

[F]I = ⟨g (z , ξ), 3[lI7(z , ξ)]⟩Y (Ω)

Up to this point, in the variational analysis, no explicit definition of X(Ω), the space used for expressing f(z,ξ), have been made. One possible representation of X7(Ω) can be obtained in terms of the two-dimensional Lagrange basis functions 37; that is,

+ ⟨fΓ (z , ξ), )[lI7(z , ξ)]⟩Y ̂(Γ)

[f]J = f J7 (z , ξ) or alternatively [f]I = f I7 (z , ξ) 16003

(58) (59)



Article

for all I, J = 0,...,7 . The indices I and J refer to the particular basis functions in the expansion of the perturbation function v(z,ξ) and the solution function f 7 (z,ξ), respectively. In an alternative interpretation, the index I refers to the collocation point, whereas the index J refers to the polynomial order of the polynomial trial function expansion 36. The Least-Squares Method as Derived through the Weighted Residual Framework (Weak Formulation). The term method of weighted residuals was originally coined by Crandall138 to describe a variety of methods (e.g., Galerkin, methods of moments, orthogonal collocation, and leastsquares) having in common the fact that minimizing the residual 9 is interpreted as equating to zero weighted integrals or inner products.139 Thus, an equivalent method for the continuous functional 1(f ; g , fΓ ) (eq 43) minimization procedure consists in finding the minimum of the inner product of the weight function w and residual function 9 , that is, ⟨w, 9⟩. Considering the generalized problem formulation 13, the minimization statement of the inner products can be given as ⟨wΩ, I , 9 Ω⟩ + ⟨wΓ, I , 9Γ⟩ = +

∫Γ wΓ,I 9Γ dΓ = 0

∫Ω wΩ,I 9 Ω dΩ + ∫Γ wΓ,I 9Γ dΓ ∫Ω

=

1 ∂ 2 ∂fI

1 ∂ 2 ∂fI

where the weight functions wI are known also as test functions. Moreover, the residuals 9 Ω and 9Γ are given as

9Γ = )f (z , ξ) − fΓ (z , ξ)

on Γ

(62)

=

∫Ω (3f (z , ξ) − g(z , ξ))2 dΩ

1 ∂ 2 ∂fI

∫Γ ()f (z , ξ) − fΓ (z , ξ))2 dΓ 7

1 ∂ 2 ∂fI

∫Ω

(3[ ∑ f J lJ7(z , ξ)] − g (z , ξ))2 dΩ J=0 7

1 ∂ 2 ∂fI

∫Γ

()[ ∑ f J lJ7(z , ξ)] − fΓ (z , ξ))2 dΓ J=0

∫Ω ( ∑ fJ 3[lJ7(z , ξ)] − g(z , ξ))3[lI7(z , ξ)]dΩ J=0

7

+

∫Γ ( ∑ fJ )[lJ7(z , ξ)] − fΓ (z , ξ)))[lI7(z , ξ)]dΓ J=0

(66)

=0

Hence, 7

∑ fJ ∫

pansion is given by eq 31, or alternatively by eq 36 in terms of global indices. The resulting residual functions are thus in terms of the solution function values f jzjξ = f J at the node points zjz and ξjξ:

Ω

J=0

3[lJ7(z , ξ)]3[lI7(z , ξ)]dΩ

7

+

∑ fJ ∫ J=0

9 Ω = 9 Ω(f0 , f1 , ..., f7 )

=

and

Γ

)[lJ7(z , ξ)])[lI7(z , ξ)]dΓ

∫Ω g(z , ξ)3[lI7(z , ξ)]dΩ + ∫Γ fΓ (z , ξ))[lI7(z , ξ)] (67)

dΓ

9Γ = 9Γ(f0 , f1 , ..., f7 )

The latter expression may alternatively be expressed by use of the bracket inner product operator:

It follows that in the least-squares method the weight functions in eq 60 are chosen as

wΓ, I =

I

7

z

wΩ, I =

∫Ω 9 Ω2 dΩ + 12 ∂∂f ∫Γ 9 2ΓdΓ

+

The spectral methods are based on using a representation of the solution f(z,ξ) throughout the domain Ω via a polynomial trial function expansion (eq 29). A nodal basis expansion in terms of the Lagrange basis polynomials is obtained by letting P P Φ Pj z (z) = l jPz(z) and Φ jξξ (ξ) = l jξξ(ξ). Thus, the series exz

I

1 ∂ = 2 ∂fI

(60)

(61)

∫Ω 9 Ω2 dΩ + 12 ∂∂f ∫Γ 9 2Γ dΓ

The latter equation may be manipulated by using the definitions of the residuals 61 and 62 and inserting the approximate solution 36. Moreover, providing that the 3 • and ) • are linear operators, the order of operation may be interchanged:

=

in Ω

∂9Γ 9 Γ dΓ ∂fI

(65)

+

9 Ω = 3f (z , ξ) − g (z , ξ)

∫Γ

=0

∫Ω wΩ,I 9 Ω dΩ

∀ I = 0, 1, ..., 7

∂9 Ω 9 Ω dΩ + ∂fI

=

∂ 9 Ω(f0 , f1 , ..., f7 ) ∂fI

∀ I = 0, 1, ..., 7

∂ 9Γ(f0 , f1 , ..., f7 ) ∂fI

∀ I = 0, 1, ..., 7

7

∑ f J ⟨3[lJ7(z , ξ)], 3[lI7(z , ξ)]⟩Y (Ω) J=0 7

(63)

+

∑ f J ⟨)[lJ7(z , ξ)], )[lI7(z , ξ)]⟩Ŷ(Γ) J=0

= ⟨g (z , ξ), 3[lI7(z , ξ)]⟩Y (Ω)

(64)

+ ⟨fΓ (z , ξ), )[lI7(z , ξ)]⟩Y ̂(Γ)

With this definition of the weight functions, the weighted residual integral statement becomes

(68)

The algebraic equation system 68 can be written on the form: 16004



Article

[L]IQ J = 3lJ7(z IQ , ξIQ )

7

∑ f J [ A ]IJ

= [ F ]I

and

[L]IQ I = 3lI7(z IQ , ξIQ ) (75)

(69)

J=0

where

[B]IQ J = )lJ7(z IQ , ξIQ )

and

[B]IQ I = )lI7(z IQ , ξIQ )

[ A ]IJ = ⟨3[lJ7(z , ξ)], 3[lI7(z , ξ)]⟩Y (Ω) + ⟨)[lJ7(z , ξ)], )[lI7(z , ξ)]⟩Y ̂(Γ)

(76)

⎡W0 ⎢ ⎢0 Λ=⎢ ⎢⋮ ⎢⎣ 0

(70)

[ F ]I = ⟨g (z , ξ), 3[lI7(z , ξ)]⟩Y (Ω) + ⟨fΓ (z , ξ), )[lI7(z , ξ)]⟩Y ̂(Γ)

(71)

Moreover, if fT = [f0 , f1 , ..., f7 ], the algebraic equation system can be expressed on the standard vector form: Af = F. This is the same results as obtained by the variational method in section. Implementation (Weak Formulation). In the following, the general implementation issues of the generalized problem formulation 3 are presented. On the basis of the variational or weighted residual principles (weak formulation), the algebraic equation system 12 (eqs 57 and 58 or alternatively eqs 70 and 71) achieved by the weak formulation of the least-squares technique can be presented as

∫Γ

+

computed by relationship 44. Alternatively, the matrix A and vector F can be written in the compact form:

method. The quadrature points and the collocation points are defined at the same locations when both type of points are determined as the roots of the same type of orthogonal

)[lJ7(z IQ , ξIQ )])[lI7(z IQ , ξIQ )]WIQ

polynomial of the same order. In this case [f]J = f J coincides

IQ = 0

with [f]IQ = f IQ.

7

=

∑

Given the boundary condition that f(z,ξ) is specified for z =

[L]IQ J [Λ]IQ IQ [L]IQ I

0, ∀ ξ, the boundary matrix B is given as

IQ = 0 7

∑

+

⎡⎡1 ⎢⎢ ⎢⎢ 0 ⎢⎢⋮ ⎢⎢ ⎢⎣ 0 ⎢ ⎢⎡ 0 ⎢⎢ ⎢⎢ 0 B = ⎢⎢⋮ ⎢⎢ ⎢⎣ 0 ⎢ ⎢ ⎢⎡ 0 ⎢⎢ ⎢⎢ 0 ⎢⎢⋮ ⎢⎢ ⎣⎢ ⎣ 0

[B]IQ J [Λ]IQ IQ [B]IQ I

IQ = 0

[F]I =

(72)

∫Ω g(z , ξ)3[lI7(z , ξ)]dΩ +

∫Γ fΓ (z , ξ))[lI7(z , ξ)]dΓ

7

=

g (z IQ , ξIQ )3[lI7(z IQ , ξIQ )]WIQ

∑ IQ = 0

7

+

∑

fΓ (z IQ , ξIQ ))[lI7(z IQ , ξIQ )]WIQ

IQ = 0 7

=

∑

(79)

approximation of the norm integrals of the least-squares

7

∑

F = LT Λg + BTΛf Γ

are normally chosen the same as the collocation points in the

)[lJ7(z , ξ)])[lI7(z , ξ)]d Γ

IQ = 0

+

(78)

To reduce time-consuming operations, the quadrature points

3[lJ7(z IQ , ξIQ )]3[lI7(z IQ , ξIQ )]WIQ

∑

A = LTΛL + BTΛB

It is noticed that a quadrature rule is applied in eqs 72 and 73.

7

=

(77)

The global quadrature weights in the diagonal matrix, Λ, is

∫Ω 3[lJ7(z , ξ)]3[lI7(z , ξ)]dΩ

[A]IJ =

... 0 ⎤ ⎥ W1 ... 0 ⎥ ⎥ ⋮ ⋱ ⋮ ⎥ 0 ... W7 ⎥⎦ 0

7

[ g ]IQ [Λ]IQ IQ [L]IQ I +

IQ = 0

∑

[ f Γ ]IQ [Λ]IQ IQ [B]IQ I

IQ = 0

(73)

[f]J = f J = f (zJ , ξJ )

0 ... 0 ⎤ ⎥ 1 ... 0 ⎥ ⋱ ⋮⎥ ⎥ 0 ... 1 ⎦ 0 ... 0 ⎤ ⎥ 0 ... 0 ⎥ ⋱ ⋮⎥ ⎥ 0 ... 0 ⎦ ⋮ 0 ... 0 ⎤ ⎥ 0 ... 0 ⎥ ⋱ ⋮⎥ ⎥ 0 ... 0 ⎦

⎡ 0 0 ... 0 ⎤ ⎡ 0 0 ... 0 ⎤ ⎤ ⎢ ⎥ ⎢ ⎥⎥ ⎢ 0 0 ... 0 ⎥ ... ⎢ 0 0 ... 0 ⎥ ⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥⎥ 0 0 ... 0 0 0 ... 0 ⎦ ⎥ ⎥ ⎡ 0 0 ... 0 ⎤ ⎡ 0 0 ... 0 ⎤ ⎥ ⎢ ⎥ ⎢ ⎥⎥ 0 0 ... 0 0 0 ... 0 ⎢ ⎥ ... ⎢ ⎥⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥ ⎢⎣ ⎥⎥ 0 0 ... 0 ⎦ 0 0 ... 0 ⎦ ⎥ ⎥ ⋮ ⋱ ⋮ ⎥ ⎡ 0 0 ... 0 ⎤ ⎡ 0 0 ... 0 ⎤ ⎥ ⎢ ⎥ ⎢ ⎥⎥ ⎢ 0 0 ... 0 ⎥ ... ⎢ 0 0 ... 0 ⎥ ⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥ ⎢⎣ ⎥⎥ 0 0 ... 0 ⎦ 0 0 ... 0 ⎦ ⎦⎥ (80)

where the submatrices are of dimension (Pξ + 1) × (Pξ + 1).

(74)

where

Moreover, the vector fΓ is given as 16005


Industrial & Engineering Chemistry Research ⎢ f (z 0 , ξ0) ⎥ ⎢ ⎥ ⎢ f (z 0 , ξ1) ⎥ ⎢ ⎥ ⎢⋮ ⎥ f Γ = ⎢ f (z , ξ )⎥ ⎢ 0 Pξ ⎥ ⎢ ⎥ ⎢0 ⎥ ⎢⋮ ⎥ ⎢⎣ ⎥⎦ 0

Article

With these definitions the weighted residual integral 89 is given as

∫Ω wΩ,I 9 Ω dΩ = ∫Ω =

1 2 ∥3f (z , ξ) − g (z , ξ)∥Y( Ω) 2

(82)

ε→ 0

d1(f (z , ξ) + εv(z , ξ); g (z , ξ)) =0 dε

∀ v(z , ξ) ∈ X(Ω)

(92)

7

(85)

[F]I = =

ξ)],

ξ)]⟩Y(Ω)

=

∑

g (z IQ , ξIQ )3[lI7(z IQ , ξIQ )]WIQ

∑

[g]IQ [Λ]IQ IQ [L]IQ I

IQ = 0

(94)

in which the matrices L and Λ are given by eqs 75 and 77. The boundary conditions can be included by manipulation of the matrix A eq 93 and vector F eq 94. Given the boundary condition that f(z, ξ) is specified for z = 0, ∀ ξ, the manipulated matrix 93 and vector 94 are given as, respectively

(87) (88)

for all I , J = 0, ..., 7 . The Least-Squares Method As Derived through the Weighted Residual Framework (Strong Formulation). Considering the generalized problem formulation 13, the inner product minimization statement is presented by eq 89, with manipulations of the resulting equation system to enforce the strong form of the boundary conditions.

∫Ω wΩ,I 9 ΩdΩ = 0

∫Ω g(z , ξ)3[lI7(z , ξ)]dΩ 7

[F]I = ⟨g (z , ξ), 3[lI7(z , ξ)]⟩Y(Ω)

⟨wΩ, I , 9 Ω⟩Ω =

(93)

IQ = 0

ξ)

3[lI7(z ,

[L]IQ J [Λ]IQ IQ [L]IQ I

7

a solvable system of algebraic equations is achieved. Following the steps of eqs 54−56, the algebraic equation system is allocated on the conventional form 12, in which ⟨3[lJ7(z ,

∑ IQ = 0

With the approximations 36 and 51 and by choosing systematically ξ), ...,

3[lJ7(z IQ , ξIQ )]3[lI7(z IQ , ξIQ )]WIQ

7

(86)

v (z , ξ ) =

∑ IQ = 0

(84)

-(g (z , ξ), fΓ (z , ξ), v(z , ξ)) = ⟨g (z , ξ), 3v(z , ξ)⟩Y(Ω)

[A]IJ =

∫Ω 3[lJ7(z , ξ)]3[lI7(z , ξ)] dΩ

=

l 77(z ,

(90)

[F]I = ⟨g (z , ξ), 3[lI7(z , ξ)]⟩Ω

=

with

l07(z ,

∀ I = 0, 1, ..., 7

(91)

[A]IJ =

We notice the similarity to eq 46 and follow the procedure of eqs 47−50. The achieved minimization statement can be presented on matrix form as

((f (z , ξ), v(z , ξ)) = ⟨3f (z , ξ), 3v(z , ξ)⟩Y(Ω)

∫Ω 9 Ω2 dΩ = 0

[A]IJ = ⟨3[lJ7(z , ξ)], 3[lI7(z , ξ)]⟩Ω

(83)

((f (z , ξ), v(z , ξ)) = -(g (z , ξ), v(z , ξ))

R Ω dΩ

Implementation (Strong Formulation). In the following, the general implementation issues of the generalized problem formulation (13) are presented. On the basis of variational or weighted residual principles (strong formulation), the algebraic equation system 12 (eqs 87 and 88) or alternatively eqs 91 and 92 achieved by the strong formulation of the least-squares technique can be presented as

with the norm defined by eq 44. The minimization statement is given as the following: Find f(z,ξ) ∈ X(Ω) such that lim

∂f I7

We notice the similarity to eq 65 and proceed in the same manner as performed in eqs 66−68. The obtained expressions for matrix A and vector F, eq 69, can thus be presented as, respectively,

(81)

The Least-Squares Method Based on Variational Principles (Strong Formulation). Considering the strong formulation of the boundary conditions, the following norm-equivalent least-squares function is defined for the generalized problem formulation (13): 1(f ; g , fΓ ) =

1 ∂ 2 ∂f I7

∂9 Ω

⎡ ⎡ 1 0 ... 0 ⎤ ⎢⎢ ⎥ ⎢ ⎢ 0 1 ... 0 ⎥ ⎢⎢⋮ ⋱ ⋮⎥ ⎢⎢ ⎥ ⎣ ⎢ 0 0 ... 1 ⎦ A′ = ⎢ [A] I = Pξ + 1, J = 0 ⎢ ⎢ ⋮ ⎢ ⎢ ⋮ ⎢ ⎢⎣ [A]I = 7, J = 0

∀ I = 0, 1, ..., 7 (89)

In the minimization statement 89, the residual 9 Ω and weighing function wΩ,I are given by eqs 61 and 63, respectively.

⎡ 0 0 ... 0 ⎤ ⎡ 0 0 ... 0 ⎤ ⎤ ⎢ ⎥ ⎢ ⎥⎥ ⎢ 0 0 ... 0 ⎥ ... ⎢ 0 0 ... 0 ⎥ ⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥⎥ 0 0 ... 0 0 0 ... 0 ⎦ ⎥ ⎥ ... ... [A]I = Pξ + 1, J = 7 ⎥ ⎥ ⎥ ⎥ ⋱ ⋮ ⎥ ... [A]I = 7, J = 7 ⎥⎦ (95)

16006


Industrial & Engineering Chemistry Research ⎡ f (z , ξ ) ⎤ ⎡[ f Γ ]I = 0 ⎤ ⎥ ⎢ 0 0 ⎥ ⎢ ⎢ f (z 0 , ξ0) ⎥ ⎢[ f Γ ]I = 1 ⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢ ⎥ ⎢ F′ = ⎢ f (z 0 , ξPξ)⎥ = ⎢[ f Γ ]I = Pξ ⎥ ⎥ ⎢ ⎥ ⎢ ⎢[F]I = Pξ + 1 ⎥ ⎢[F]I = Pξ + 1⎥ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢⋮ ⎥ ⎢⋮ ⎥ ⎢[F] ⎥ ⎢ ⎣ I = P ⎦ ⎣[F]I = P ⎦

Article

fist column and rows are removed from matrix A 95, and the Pξ + 1 fist rows are removed from the vector F 96. The vector F is further manipulated with the A matrix coefficients to enforce the boundary conditions:

(96)

It should be noticed that the global numbering procedure as sketched in Figure 15 is adopted. Galerkin and Tau Methods. Consider the generalized problem formulation defined by eqs 13a and 13b. The inner product minimization statement is given by eq 89. In both the Galerkin and tau methods, the weighting functions are chosen identical to the basis function, that is, wΩ, I = lI7(z , ξ) 7 ⟨lI7(z , ξ), 9 Ω (z , ξ , f 07 , .., f J7 , .., f 77 )⟩

∫Ω 9 Ω7(z , ξ , f 07 , .., f J7 , .., f 77 )lI7(z , ξ)dΩ 7

=

∫Ω

[ ∑ f J7 3lJ7(z , ξ) − g (z , ξ)]lI7(z , ξ)dΩ J=0

7

=

7

J=0

7

=

J=0

(101)

=

7

∑ f J7 ∑

⎡[F] ⎤ ⎢ I = Pξ + 1 − f (z 0 , ξ0)[A]I = Pξ + 1, J = 0 ... ⎥ ⎢ − f (z 0 , ξP )[A]I = P + 1, J = P + 1 ⎥ ξ ξ ξ ⎢ ⎥ ⎥ F″ = ⎢ ⋮ ⎢ ⎥ ⎢[F]I = 7 − f (z 0 , ξ0)[A]I = 7, J = 0 ... ⎥ ⎢ ⎥ ⎢⎣ − f (z 0 , ξPξ)[A]I = 7, J = Pξ + 1 ⎥⎦

7 (z , ξ , f 07 , .., f J7 , .., f 77 )⟩ ⟨δ(x − xI ), 9 Ω

∑ [ ∑ f J7 3lJ7(z IQ , ξIQ ) − g(z IQ , ξIQ )]lI7(z IQ , ξIQ )WIQ IQ = 0

(100)

Orthogonal Collocation Method. Consider the generalized problem formulation defined by eqs 13a and 13b. The inner product minimization statement is given by eq 89. In the orthogonal collocation method, the weighting function is the Dirac delta function, that is, ωΩ,I(x) = δ(x − xI). Accordingly, at the location x = xI the residual is forced to zero such that the generalized problem (13) is satisfied at these points. Thus the inner products are given as

Hence,

=

⎡[A]I = P + 1, J = 0 ... ... [A]I = P + 1, J = 7 ⎤ ξ ξ ⎢ ⎥ ⎢ ⎥ ⋮ ⎥ A″ = ⎢ ⎢ ⎥ ⋮ ⋱ ⋮ ⎢ ⎥ ⎢ [A] ⎥ ... [ A ] I = 7, J = 0 I = 7, J = 7 ⎦ ⎣

lI7(z IQ , ξIQ )3lJ7(z IQ , ξIQ )WIQ

7 (z , ξ , f 07 , ..., f J7 , ..., f P7 )|x = xI = 9Ω

IQ = 0

7

−

∫Ω 9 Ω7(z , ξ , f 07 , ..., f J7 , ..., f P7 )δI dΩ 7

lI7(z IQ , ξIQ )g (z IQ , ξIQ )WIQ = 0

∑

∀ I = 0, 1, ..., 7

=

IQ = 0

∑ f J7 3lJ7(zI , ξI ) − g(zI , ξI ) = 0

∀ I = 0, 1, ..., 7

J=0

(97)

(102)

in which 9 Ω is defined by eq 61. Moreover, the polynomial trial function expansion 29 has been adopted to approximated the solution function. The latter equation can be presented on the conventional matrix form 12:

where the polynomial trial function expansion (31) has been applied to approximate the solution function. The approximated solution f 7 (zI , ξI ) for I = 0, 1, ..., 7 obtained by the collocation method is thus represented by finding a particular set of unknown coefficients {f J7 , J = 0, 1, ..., P} which satisfy

7

[A]IJ =

lI7(z IQ , ξIQ )3lJ7(z IQ , ξIQ )WIQ

∑

eqs 102. Equation 102 can be presented on the form 12:

IQ = 0 7

=

lI7(z IQ ,

∑

ξIQ )[L]IQ J WIQ

IQ = 0

(98)

7

[F]I =

∑ IQ = 0

lI7(z IQ , ξIQ )g (z IQ , ξIQ )WIQ

[ A ]IJ = 3lJ7(zI , ξI ) = [L]IJ

(103)

[ F ]I = g (zI , ξI )

(104)

f = [ f 07 , f 17 , ..., f 77 ]T

(105)

in which L is defined in accordance to eq 75. Further, the treatment of the strong form implementation of the boundary conditions can follow the approach of eqs 100 and 101. However, due to the simplicity, the approach in eqs 95 and 96 is employed in the present study for the orthogonal collocation method. Implementation Issues for Spectral Solution of the PBE. Further details on the implementation issues using the leastsquares method for the solution of the PBE 2 is given by Solsvik and Jakobsen.98 The implementation framework

(99)

and with fT = [f 07 , f 17 , ..., f 77 ]T . Moreover, L is defined in accordance to eq 75. The equation system must be manipulated in order to enforce the boundary conditions. Given the boundary condition that f(z,ξ) is specified for z = 0, ∀ ξ. For the tau method, the strong form implementation of the boundary conditions are handled in the same manner as in eqs 95 and 96. On the other hand, the Galerkin method requires relative more comprehensive manipulations. Thus, the Pξ + 1 16007



Article

discussed by Solsvik and Jakobsen98 can be adopted for the Galerkin, tau, and orthogonal collocation methods noting the differences in the algebra of the solution techniques, as presented in the Appendix. The problem operator matrix L and source vector g that are applied to obtain the equation system 12 of the orthogonal collocation (eqs 103 and 104), Galerkin (eqs 98 and 99), tau (eqs 98 and 99), and leastsquares (strong form) (eqs 93 and 94) methods are given in the sequence. Manipulation of the resulting equation system to enforce the boundary conditions follows the technique as presented by eqs 95 and 96 (orthogonal collocation, tau, strong form least-squares) and 102 and 103 (Galerkin). In addition to the problem operator matrix L and source vector g, the boundary operator matrix B and boundary values fΓ constitute the weak formulation of the least-squares method (eqs 72 and 73). Boundary Conditions. Considering the PB problem 2, a set of boundary conditions are required. In the present study, the following set of boundary conditions are employed: (i) f(z,ξ) = 0 for ξ = ξmin, ∀ z, and (ii) f(z,ξ) = specified for z = 0, ∀ ξ. The specified distribution function at z = 0 is given in, for example, Figure 12. The Linear System of Algebraic Equations. To obtain a system of algebraic equations, the PBE 2 and the boundary conditions must be transformed into a discrete form. In spectral methods, the solution function is approximated in terms of a polynomial trial function expansion:

[ L1]IQ J = vz

[ L 2 ]IQ J = −

×

jz = 0 jξ = 0

∂f (z , ξ) ≈ ∂z

∂f (z , ξ) ≈ ∂ξ

∑ ∑ f j7(J)j (J) jz = 0 jξ = 0

∑∑ jz = 0 jξ = 0

ξiqξ(IQ )vz

∂ρd (z)

3ρd (z iqz(IQ ))

∂z

iqz(IQ )

∂ Pξ l (ξ ) ∂ξ jξ (J )

l jPz(J )(z iqz(IQ )) z

(111)

P

z

ξ

∂ρd (z) vz 3ρd (z iqz(IQ )) ∂z

iqz(IQ )

(112)

Breakage death:

ξ

ξ

∂l jPz(J )(z) z

∂z

[ L4 ]IQ J = H(ξiqξ(IQ ) −

3

3 2ξmin ) P

× b(ξiqξ(IQ ))l jPz(J )(z iqz(IQ ))l j ξ(J )(ξiqξ(IQ )) z

(113)

ξ

P

l j ξ(J )(ξ) ξ

Breakage birth:

(107) P

Pξ

Pz

(110)

(106)

Pξ z

ξ

[ L3 ]IQ J = −l j ξ(J )(ξiqξIQ )l jPz(J )(z iqz(IQ ))

The differentiation of the discrete solution approximation may be presented as Pz

l j ξ(J )(ξiqξ(IQ ))

∂z

iqξ(IQ )

ξ

z

z

ξ

P

z

Convection in the property coordinate ξ, i.e. bubble growth:

∑ ∑ f j7(J)j (J) l jP(J)(z)l jP(J)(ξ) z

∂l jPz(J )(z) iqz(IQ )

Pξ

Pz

f (z , ξ ) ≈

Convection in the physical space z:

f j7(J )j (J ) l jPz(J )(z) z z ξ

∂l j ξ(J )(ξ)

3 3 [ L5 ]IQ J = −H( 3 ξmax − ξmin − ξiqξ(IQ ))

ξ

∂ξ

× V (ξiqξ(IQ ))l jPz(J )(z iqz(IQ ))

(108)

z

P ⎡ ⎤ l j ξ(J )(ζiqζ ) ⎢ ⎥ Pζ ξ (W iqζ ξiqξ(IQ ) , ξmax )⎥ × ∑ ⎢h(ξiqξIQ , ζiqζ )b(ζiqζ ) V (ζiqζ ) iqζ = 0 ⎢ ⎥⎦ ⎣

The PBE 2 is an integro−differential equation. Thus, appropriate quadrature rules are required for the numerical solution. A problem arises as the solution vector f is only available in the collocation points, f jzjξ, and not necessarily in the ζ-points, f jzjζ (see the integral terms for breakage birth and coalescence birth and death in eq 2). Hence,

∫ξ

ξ2

Pζ

(114)

Coalescence death:

f (z , ζ )d ζ

1

Piqζ

≈

Pz

3 3 [ L6 ]IQ J = H( 3 ξmax − ξmin − ξiqξ(IQ ))f jν(−J )1j (J )

Pξ

∑∑∑ iqζ = 0 jz = 0 jξ = 0

P P f j7(J )j (J ) l jPz(J )(z)l j ξ(J )(ζiqζ )W iqζζ(ξ1 , z ξ z ξ

z

ξ2) ×

(109)

Considering the PBE 2, the problem matrix [L], as defined by eq 75, is given by eqs 110−117.

l jPz(J )(z iqz(IQ )) z ρd (z iqz(IQ ))

ξ

P

Pζ

×

∑

[c(ξiqξ(IQ ) , ζiqζ )

l j ξ(ζi )(ζiqζ ) ξ

iqζ = 0

qζ

V (ζiqζ )

P

× W iqζζ(ξmin , [(ξmax )3 − (ξiqξ(IQ ))3 ]1/3 )] 16008

(115)



Article

⎡⎡1 ⎢⎢ ⎢⎢ 0 ⎢⎢⋮ ⎢⎢ ⎢⎣ 0 ⎢ ⎢⎡ 0 ⎢⎢ 0 ⎢⎢ B = ⎢⎢⋮ ⎢ ⎢⎣ 0 ⎢ ⎢ ⎢⎡ ⎢⎢ 0 ⎢⎢ 0 ⎢⎢ ⎢⎢⋮ ⎣⎢ ⎣ 0

Coalescence birth: [ L 7 ]IQ J = −H(ξiqξ(IQ ) − ×

3

3 2ξmin )

P ξiq2ξ(IQ )V (ξiqξ(IQ )) l j z(J )(z iqz(IQ )) z

ρd (z iqz(IQ ))

2

⎡ 3 3 1/3 ⎢ c([[ξiqξ(IQ )] − [ζiqζ ] ] , ζiqζ ) × ∑⎢ [[ξiqξ(IQ )]3 − [ζiqζ ]3 ]2/3 iqζ ⎣ ⎢ ×

×

Pξ 1 3 3 1/3 f jν− (J )j l j (J )([[ξiqξ(IQ )] − [ζiqζ ] ] ) z

ζ

ξ

V (ζiqζ ) P W iqζζ(ξmin ,

[V (ξiqξ(IQ )) − V (ζiqζ )] 3

3 1/3

[[ξiqξ(IQ )] − [ξmin] ]

⎤ ⎥ )⎥ ⎥⎦

⋮ 0 ... 0 ⎤ ⎥ 0 ... 0 ⎥ ⋱ ⋮⎥ ⎥ 0 ... 0 ⎦

where the submatrices are of dimension (Pξ + 1) × (Pξ + 1). Moreover, the boundary condition vector is given by [ f Γ ]IQ = [f (z 0 , ξ0), f (z 0 , ξ1), ..., f (z 0 , ξPξ), f (z1 , ξ0), 0, ...,

7

■

∑ [ Lq ]I J Q

(117)

q=1

0⎤ ⎥ 0⎥ ⋱ ⋮⎥ ⎥ 0 0⎦ 0 0

⎡ 0 0 ... 0 ⎤ ⎡ 0 0 ... 0 ⎤ ⎤ ⎢ ⎥ ⎢ ⎥⎥ ⎢ 0 0 ... 0 ⎥ ... ⎢ 0 0 ... 0 ⎥ ⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥⎥ 0 0 ... 0 0 0 ... 0 ⎦ ⎥ ⎥ ⎡1 0 ⎡ 0 0 ... 0 ⎤ ⎥ 0⎤ ⎢ ⎥ ⎢ ⎥ 0⎥ 0 0 ... 0 ⎥ ⎥ ⎢0 0 ... ⎢ ⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥ ⎢⎣ ⎥ 0 0 0⎦ 0 0 ... 0 ⎦ ⎥ ⎥ ⋮ ⋱ ⋮ ⎥ ⎡ 0 0 ... 0 ⎤ ⎡ 1 0 ... 0 ⎤ ⎥⎥ ⎢ ⎥ ⎢ ⎥ ⎢ 0 0 ... 0 ⎥ ... ⎢ 0 0 ... 0 ⎥ ⎥⎥ ⎢⋮ ⎢⋮ ⋱ ⋮⎥ ⋱ ⋮⎥⎥ ⎢⎣ ⎥⎦ ⎢⎣ ⎥ 0 0 ... 0 0 0 ... 0 ⎦ ⎦⎥ (121)

(116)

The total problem matrix is obtained by summarizing the submatrices:

[L]IQ J =

0 ... 0 ⎤ ⎥ 1 ... 0 ⎥ ⋱ ⋮⎥ ⎥ 0 ... 1 ⎦

f (z 2 , ξ0), 0, ..., ..., f (z Pz , ξ0), 0, ...)]T

(122)

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]; hugo.jakobsen@ chemeng.ntnu.no.

In the above equations, the heaviside function is defined with the following properties:

Notes

The authors declare no competing financial interest.

⎧ = 0 for Ψ < 0 H(Ψ)⎨ ⎩ = 1 for Ψ ≥ 0

■ ■

(118)

ACKNOWLEDGMENTS The Research Council of Norway (NFR) is thanked for financial support (J. Solsvik).

Thus, with the heaviside function special care is introduced at the boundaries of the ξ-domain in order to preserve the mass of the system.

Latin Letters

In the present study, the Picard iteration technique is

b = Breakage frequency (1/s) c = Coalescence closure f = Solution function; density probability function fΓ = Solution function at the boundary of the Ω-domain fd,m (ξ, r, t) = Length based mass density function (kg/m3m) i = Local index associated with the collocation points I = Global index associated with the collocation points j = Local index associated with the polynomial order in the truncated series expansion J = Global index associated with the polynomial order in the truncated series expansion J = Jacobian Jα,β j = Jacobi polynomial of degree j g = Gravitational acceleration (m/s2) H (ψ) = Heaviside function h(ξ, ζ) = Daughter size redistribution function (1/m) ho = Initial film thickness (m) hf = Critical film thickness for rapture (m) KC = Prefactor k1 = Breakage kernel parameter (s) k2 = Breakage kernel parameter (−) L = Column height (m) N = Legendre polynomial

employed in the nonlinear coalescence birth and death terms. 1 The term f j7(,Jν− in the coalescence birth term is approximated )j z

ζ

by Pξ

f jν(−J )1j = z

ζ

∑ f jν(−J)1j jξ

z

ξ

P

× l jξξ(ξiqζ ) (119)

The source term vector is given as

[ g ]IQ = 0

NOTATION

(120)

Considering the weak formulation of the least-squares method. On the basis of the boundary conditions presented, the boundary condition matrix B, defined by eq 76, is given as 16009



Article

ξ = Property coordinate ξ Y (Ω) = Space for the norms measuring the residual in Ω Ŷ (Γ) = Space for the norms measuring the residual in Γ

Pz = Polynomial order in z coordinate Pqz = Order of quadrature Pξ = Polynomial order in the property coordinate ξ 7 = Order of the multidimensional problem p = Pressure (Pa) p0 = Atmospheric pressure (Pa) re = Equivalent radius (m) r = Coordinate in physical space (m) 9 = Residual :BD(ξ , r, t ) = Source term due to breakage and coalescence V = Volume (m3) vr = Velocity in physical space (m/s) υξ = Growth velocity (m/s) υz = Velocity in physical space, z-direction (m/s) υ = Perturbation function W = Quadrature weight w = Weighting function X (Ω) = Trial solution space Y (Ω) = Space for the norms measuring the residual in Ω Ŷ (Γ) = Space for the norms measuring the residual on Γ z = Spatial coordinate (m)

Superscript

d = Dimension ν = iteration 7 = Polynomial order Pz = Polynomial order in z-coordinate Pξ = Polynomial order in ξ-coordinate T = Transpose

■

REFERENCES

(1) Deckwer, W. D. Bubble Column Reactors; John Wiley & Sons: New York, 1992. (2) Jakobsen, H. A. Chemical Reactor Modeling: Multiphase Reactive Flows; Springer: Berlin, 2008. (3) Jakobsen, H. A.; Sannæs, B. H.; Grevskott, S.; Svendsen, H. F. Modeling of vertical bubble-driven flows. Ind. Eng. Chem. Res. 1997, 36, 4052−4074. (4) Jakobsen, H. A.; Lindborg, H.; Dorao, C. A. Modeling of bubble column reactors: Progress and limitations. Ind. Eng. Chem. Res. 2005, 44, 5107−5151. (5) Dorao, C. A. High order methods for the solution of the population balance equation with applications to bubbly flows. Ph.D. Thesis, Norwegian University of Science and Technology (NTNU), 2006. (6) Ramkrishna, D. Population balance: Theory and Applications to Particulate Systems in Engineering; Academic Press: San Diego, 2000. (7) Ramkrishna, D.; Mahoney, A. W. Population balance modeling. Promise for the future. Chem. Eng. Sci. 2002, 57, 595−606. (8) Ramkrishna, D. The status of population balances. Rev Chem Eng 1985, 3, 49−95. (9) Randolph, A. D.; Larson, M. A. Theory of particulate processes. Analysis and techniques of continuous crystallization, 2nd ed.; Academic Press: San Diego, 1988. (10) Sporleder, F.; Borka, Z.; Solsvik, J.; Jakobsen, H. A. On the population balance equation. Rev. Chem. Eng. 2012, 28, 149−169. (11) Rigopoulos, S. Population balance modelling of polydispersed particles in reactive flows. Prog. Energy Combust. Sci. 2010, 36, 412− 443. (12) Ramkrishna, D. Analysis of population balance-IV: The precise connection between Monte Carlo simulation and population balances. Chem. Eng. Sci. 1981, 36, 1203−1209. (13) Mahoney, A. W.; Ramkrishna, D. Efficient solution of population balance equations with discontinuities by finite elements. Chem. Eng. Sci. 2002, 57, 1107−1119. (14) Dorao, C. A.; Jakobsen, H. A. The quadrature method of moments and its relationship with the method of weighted residuals. Chem. Eng. Sci. 2006, 61, 7795−7804. (15) Singh, P. N.; Ramkrishna, D. Solution of population balance equations by MWR. Comput. Chem. Eng. 1977, 1, 23−31. (16) Dorao, C. A.; Jakobsen, H. A. A least squares method for the solution of population balance problems. Comput. Chem. Eng. 2006, 30, 535−547. (17) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-I. A fixed pivot technique. Chem. Eng. Sci. 1996, 51, 1311−1332. (18) Marchisio, D. L.; Fox, R. O. Solution of population balance equations using the direct quadrature method of moments. J. Aerosol Sci. 2005, 36, 43−73. (19) Marchisio, D. L.; Pikturna, J. T.; Fox, R. O.; Vivil, R. D.; Barresi, A. A. Quadrature method of moments for population-balance equations. AIChE 2003, 49, 1266−1276. (20) Marchisio, D. L.; Vigil, R. D.; Fox, R. O. Quadrature method of moments for aggregation-breakage processes. J. Colloid Interface Sci. 2003, 258, 322−334.

Greek Letters

α = Basis coefficients in truncated series expansion α̅ i = Mean volume fraction of liquid Φ = Basis (or trial) functions in truncated series expansion ξ = Bubble diameter (m) ζ = Bubble diameter (m) ρ = Density (kg/m3) l = Lagrange polynomial Ω = Domain 3 = Linear operator 1 = Functional ) = Linear boundary operator Γ = Boundary σ = Gas−liquid surface tension (N/m) ε = Energy dissipation (m2/s3) ε = Amplitude variable in variational analysis β = Parameter in coalescence kernel Λ = Quadrature weight matrix ω = Weight function in the weighted residual principle Subscript

B = Birth b = Breakage C = Coalescence c = Continuous phase D = Death d = Dispersed phase iz = Local index in z-coordinate iξ = Local index in ξ-coordinate iqξ(IQ) = Local index in quadrature point ξ iqζ = Local index in quadrature point ζ iqz(IQ) = Local index in quadrature point z I = Global index IQ = Global index jz = Local index in z-coordinate jξ = Local index in ξ-coordinate J = Global index associated with the polynomials of the truncated series expansion m = Mass basis min = Minimum max = Maximum z = Physical space coordinate z 16010



Article

(42) Han, L.; Gong, S.; Li, Y.; Ai, Q.; Luo, H.; Liu, Z.; Liu, Y. A novel theoretical model of breakage rate and daughter size distribution for droplet in turbulent flows. Chem. Eng. Sci. 2013, 102, 186−199. (43) Raikar, N. B.; Bhatia, S. R.; Malone, M. F.; Henson, M. A. Experimental studies and population balance equation models for breakage prediction of emulsion drop size distributions. Chem. Eng. Sci. 2009, 64, 2433−2447. (44) Mallikarjunan, V.; Pushpavanam, S.; Immanuel, C. D. Parameter estimation strategies in batch emulsion polymerization. Chem. Eng. Sci. 2010, 65, 4967−4982. (45) Jildeh, H. B.; Attarakih, M.; Mickler, M.; Bart, H. J. Can. J. Chem. Eng., DOI: 10.1002/cjce.21892. (46) Singh, K. K.; Mahajani, S. M.; Shenoy, K. T.; Ghosh, S. K. Population balance modeling of liquid-liquid dispersions in homogeneous continuous-flow stirred tank. Ind. Eng. Chem. Res. 2009, 48, 8121−8133. (47) Laakkonen, M.; Alopaeus, V.; Aittamaa, J. Validation of bubble breakage, coalescence and mass transfer models for gas-liquid dispersion in agitated vessel. Chem. Eng. Sci. 2006, 61, 218−228. (48) Laakkonen, M.; Moilanen, P.; Alopaeus, V.; Aittamaa, J. Modelling local bubble size distributions in agitated vessels. Chem. Eng. Sci. 2007, 62, 721−740. (49) Becker, P. J.; Puel, F.; Henry, R.; Sheibat-Othman, N. Investigation of discrete population balance models and breakage kernels for dilute emulsification systems. Ind. Eng. Chem. Res. 2011, 50, 11358−11374. (50) Becker, P. J.; Puel, F.; Chevalier, Y.; Sheibat-Othman, N. Can. J. Chem. Eng., DOI: 10.1002/cjce.21885. (51) Lucas, D.; Krepper, E.; Prasser, H. M. Development of cocurrent air−water flow in a vertical pipe. Int. J. Multiphase Flow 2005, 31, 1304−1328. (52) Prasser, H. M.; Beyer, M.; Carl, H.; Gregor, S.; Lucas, D.; Pietruske, H.; Schutz, P.; Weiss, F. P. Evolution of the structure of a gas-liquid two-phase flow in a large vertical pipe. Nucl. Eng. Des. 2007, 237, 1848−1861. (53) Lucas, D.; Beyer, M.; Kussin, J.; Schütz, P. Benchmark database on the evolution of two-phase flows in a vertical pipe. 2008; XCFD4NRS, Experiments and CFD code applications to nuclear reactor safety. (54) Duan, X. Y.; Cheung, S. C. P.; Yeoh, G. H.; Tu, J. Y.; Krepper, E.; Lucas, D. Gas-liquid flows in medium and large vertical pipes. Chem. Eng. Sci. 2011, 66, 872−883. (55) Krepper, E.; Lucas, D. Population balance model for the CFD simulation of adiabatic and diabatic two phase gas liquid flows. 2012. Unpublished work. (56) Dorao, C. A.; Lucas, D.; Jakobsen, H. A. Prediction of the evolution of the dispersed phase in bubbly flow problems. Appl. Math. Model. 2008, 32, 1813−1833. (57) Zhu, Z.; Dorao, C. A.; Lucas, D.; Jakobsen, H. A. On the coupled solution of a combined population balance model using the least-squares spectral element method. Ind. Eng. Chem. Res. 2009, 48, 7994−8006. (58) Maaß, S.; Gäbler, A.; Zaccone, A.; Paschedag, A. R.; Kraume, M. Experimental investigations and modelling of breakage phenomena in stirred liquid/liquid systems. Chem. Eng. Res. Des. 2007, 85, 703−709. (59) Zaccone, A.; Gäbler, A.; Maaß, S.; Marchisio, D.; Kraume, M. Drop breakage in liquid−liquid stirred dispersions: Modeling of single drop breakage. Chem. Eng. Sci. 2007, 62, 6297−6307. (60) Maaß, S.; Kraume, M. Determination of breakage rates using single drop experiments. Chem. Eng. Sci. 2012, 70, 146−164. (61) Buffo, A.; Vanni, M.; Marchisio, D. L.; Fox, R. O. Multivariate quadrature-based moments methods for turbulent polydisperse gas− liquid systems. Int. J. Multiphase Flow 2013, 50, 41−57. (62) Buffo, A.; Vanni, M.; Marchisio, D.; Fox, R. O. Comparison between different methods for turbulent gas-liquid systems by using multivariate population balances. 2012; 8th International Conference on CFD in Oil & Gas, Metallurgical and Process Industries, SINTEF/ NTNU, Trondheim Norway.

(21) Tambour, Y.; Seinfeld, J. H. Sectional representations for simulating aerosol dynamics. J. Colloid Interface Sci. 1980, 76, 541− 556. (22) Yeoh, G. H.; Tu, J. Y. Two-fluid and population balance models for subcooled boiling flow. Appl. Math. Model. 2006, 30, 1370−1391. (23) Dorao, C. A.; Fernandino, M.; Patruno, L. E.; Dupuy, P. M.; Jakobsen, H. A.; Svendsen, H. F. Macroscopic description of dropletfilm interaction for gas-liquid systems. Appl. Math. Model. 2009, 33, 3309−3318. (24) Vafa, E.; Shahrokhi, M.; Abedini, H. Solution of population balance equations in emulsion polymerization using method of moments. Chem. Eng. Commun. 2013, 200, 20−49. (25) Balakin, B. V.; Hoffman, A. C.; Kosinski, P. Population balance model for nucleation, growth, aggregation, and breakage of hydrate particles in turbulent flow. AIChE J. 2010, 56, 2052−2062. (26) Gerstlauer, A.; Mitrović, A.; Motz, S.; Gilles, E. D. A population balance model for crystallization process using two independent particle properties. Chem. Eng. Sci. 2001, 56, 2553−2565. (27) Poon, J. M. H.; Immanuel, C. D.; Doyle, F. J.; Lister, J. D. A three-dimensional population balance model of granulation with a mechanistic representation of the nucleation and aggregation phenomena. Chem. Eng. Sci. 2008, 63, 1215−1329. (28) Celnik, M.; Patterson, R.; Kraft, M.; Wagner, W. Coupling a stochastic soot population balance to gas-phase chemistry using operator splitting. Combust. Flame 2007, 148, 158−176. (29) Fredrickson, A. G.; Mantzaris, N. V. A new set of population balance equations for microbial and cell cultures. Chem. Eng. Sci. 2002, 57, 2265−2278. (30) Solsvik, J.; Tangen, S.; Jakobsen, H. A. On the constitutive equations for fluid particle breakage. Rev. Chem. Eng. 2013, 29, 241− 356. (31) Silva, L. F. L. R.; Rodrigues, R. C.; Mitre, J. F.; Lage, P. L. C. Comparison of the accuracy and performance of quadrature-based methods for population balance problems with simultaneous breakage and aggregation. Comput. Chem. Eng. 2010, 34, 286−297. (32) JunWei, S.; ZhwoLin, G.; Yun, X. X. Advances in numerical methods for the solution of population balance equations for dispersed phase systems. Sci. China Ser. B, Chem. 2009, 52, 1063−1079. (33) Mantzaris, N. V.; Daoutidis, P.; Srienc, F. Numerical solution of multi-variable cell population balance models. II. Spectral methods. Comput. Chem. Eng. 2001, 25, 1441−1462. (34) Attarakih, M.; Al-Zyod, S.; Abu-Khader, M.; Bart, H. J. PPBLAB: A new multivariate population balance environment for particulate system modelling and simulation. Proc. Eng. 2012, 42, 1445−1462. (35) Braumann, A.; Man, P. L. W.; Kraft, M. Statistical approximation of the inverse problem in multivariate population balance modeling. Ind. Eng. Chem. Res. 2010, 49, 428−438. (36) Ramachandran, R.; Barton, P. I. Effective parameter estimation within a multi-dimensional population balance model framework. Chem. Eng. Sci. 2010, 65, 4884−4893. (37) Sajjadi, B.; Raman, A. A. A.; Shah, R. S. S. R. E.; Ibrahim, S. Review on applicable breakup/coalescence models in turbulent liquidliquid flows. Rev. Chem. Eng. 2013, 29, 131−158. (38) Borka, Z.; Jakobsen, H. A. Evaluation of breakage and coalescence kernels for vertical bubbly flows using a combined multifluid-population balance model solved by least squares method. Proc. Eng. 2012, 42, 623−633. (39) Borka, Z.; Jakobsen, H. A. On the modeling and simulation of higher order breakage for vertical bubbly flows using the least squares method: Application for bubble column and pipe flows. Proc. Eng. 2012, 42, 1270−1281. (40) Solsvik, J.; Borka, Z.; Becker, P. J.; Sheibat-Othman, N.; Jakobsen, H. A. Can. J. Chem. Eng., DOI: 10.1002/cjce.21876. (41) Nguyen, V. T.; Song, C. H.; Bae, B. U.; Euh, D. J. Modeling of bubble coalescence and break-up considering turbulent suppression phenomena in bubbly two-phase flow. Int. J. Multiphase Flow 2013, 54, 21−42. 16011



Article

(63) Bhole, M. R.; Joshi, J. B.; Ramkrishna, D. CFD simulation of bubble columns incorporating population balance modeling. Chem. Eng. Sci. 2008, 63, 2267−2282. (64) Hagesaether, L.; Jakobsen, H. A.; Svendsen, H. F. Modeling of the dispersed-phase size distribution in bubble columns. Ind. Eng. Chem. Res. 2002, 41, 2560−2570. (65) Sha, Z.; Laari, A.; Turunen, I. Multi-phase-multi-size-group model for the inclusion of population balances into the CFD simulation of gas-liquid bubbly flows. Chem. Eng. Technol. 2006, 29, 550−559. (66) Chen, P.; Dudukovic̀, M. P.; Sanyal, J. Three-dimensional simulation of bubble column flows with bubble coalescence and breakup. AIChE J. 2005, 51, 696−712. (67) Chen, P.; Sanyal, J.; Dudukovic̀, M. P. CFD modeling of bubble columns flows: implementation of population balance. Chem. Eng. Sci. 2004, 59, 5201−5207. (68) Krepper, E.; Lucas, D.; Frank, T.; Prasser, H.-M.; Zwart, P. J. The inhomogeneous MUSIG model for the simulation of polydispersed flows. Nucl. Eng. Des. 2008, 238, 1690−1702. (69) Sanyal, J.; Marchisio, D. L.; Fox, R. O.; Dhanasekharan, K. On the comparison between population balance models for CFD simulations of bubble columns. Ind. Eng. Chem. Res. 2005, 44, 5063−5072. (70) Morud, J. Implementation of the sectional quadrature method of moments in Fluent. 2011. (71) Zhu, Z. The least-squares spectral element method solution of the gas-liquid multi-fluid model coupled with the population balance equation. Ph.D. Thesis, Norwegian University of Science and Technology (NTNU), 2009. (72) Patruno, L. E. Experimental and numerical investigations of liquid fragmentation and droplet generation for gas processing at high pressures. Ph.D. Thesis, Norwegian University of Science and Technology (NTNU), 2010. (73) Sporleder, F. Simulation of chemical reactors using the leastsquares spectral element method. Ph.D. Thesis, Norwegian University of Science and Technology (NTNU), 2011. (74) Nayak, A. K.; Borka, Z.; Patruno, L. E.; Sporleder, F.; Dorao, C. A.; Jakobsen, H. A. A combined multifluid-population balance model for vertical gas-liquid bubble-driven flows considering bubble column operating conditions. Ind. Eng. Chem. Res. 2011, 50, 1786−1798. (75) Dorao, C. A.; Jakobsen, H. A. Numerical calculation of the moments of the population balance equation. J. Comput. Appl. Math. 2006, 196, 619−633. (76) Lafi, A. Y.; Reyes, J. N. General particle transport equation. Final report 1994. (77) Lasheras, J. C.; Eastwood, C.; Martinez-Bazan, C.; Montanes, J. L. A review of statistical models for the break-up of an immisible fluid immersed into a fully developed turbulent flow. Int. J. Multiphase Flow 2002, 28, 247−278. (78) Lathouwers, D.; Bellan, J. Yield optimization and scaling of fluidized beds for tar production from biomass. Energy Fuels 2001, 15, 1247−1262. (79) Dorao, C. A.; Jakobsen, H. A. Least-squares spectral method for solving advective population balance problems. J. Comput. Appl. Math. 2007, 201, 247−257. (80) Dorao, C. A.; Jakobsen, H. A. Time-property least-squares spectral method for population balance equations. J. Math. Chem. 2009, 46, 770−780. (81) Dorao, C. A.; Jakobsen, H. A. hp-adaptive least squares spectral element method for population balance equations. Appl. Num. Math. 2008, 58, 563−576. (82) Dorao, C. A.; Jakobsen, H. A. Time-space-property least squares spectral method for population balance problems. Chem. Eng. Sci. 2007, 62, 1323−1333. (83) Patruno, L. E.; Dorao, C. A.; F, S. H.; Jakobsen, H. A. Analysis of breakage kernels for population balance modelling. Chem. Eng. Sci. 2009, 64, 501−508.

(84) Zhu, Z.; Dorao, C. A.; Jakobsen, H. A. Solution of bubble number density with breakage and coalescence in a bubble column by least-squares method. Prog. Comput. Fluid Dyn. 2009, 9, 436−446. (85) Zhu, Z.; Dorao, C. A.; Jakobsen, H. A. A least-squares method with direct minimization for the solution of the breakage-coalescence population balance equation. Math. Comput. Simul. 2008, 79, 716− 727. (86) Zhu, Z.; Dorao, C. A.; Jakobsen, H. A. Mass conservative solution of the population balance equation using the least-squares spectral element method. Ind. Eng. Chem. Res. 2010, 49, 6204−6214. (87) Sporleder, F.; Dorao, C. A.; Jakobsen, H. A. Model based on population balance for the simulation of bubble columns using methods of the least-squares type. Chem. Eng. Sci. 2011, 66, 3133− 3144. (88) Borka, Z.; Jakobsen, H. A. Least squares higher order method for the solution of a combined multifluid-population balance model: Modeling and implementation issues. Proc. Eng. 2012, 42, 1121−1132. (89) Solsvik, J.; Jakobsen, H. A. A numerical study of a two property catalyst/sorbent pellet design for the sorption-enhanced steammethane reforming process: Modeling complexity and parameter sensitivity study. Chem. Eng. J. 2011, 178, 407−422. (90) Solsvik, J.; Jakobsen, H. A. Modeling of multicomponent mass diffusion in porous spherical pellets: Application to steam methane reforming and methanol synthesis. Chem. Eng. Sci. 2011, 66, 1986− 2000. (91) Solsvik, J.; Tangen, S.; Jakobsen, H. A. On the consistent modeling of porous catalyst pellets: Mass and molar formulations. Ind. Eng. Chem. Res. 2012, 51, 8222−8236. (92) Van den Bosch, B.; Padmanabhan, L. Use of orthogonal collocation methods for the modeling of catalyst particles-I. Analysis of the multiplicity of the solutions. Chem. Eng. Sci. 1974, 29, 1217−1225. (93) Lefévre, L.; Dochain, D.; Feyo de Azevedo, S.; Magnus, A. Optimal selection of orthogonal polynomials applied to the integration of chemical reactor equations by collocation methods. Comput. Chem. Eng. 2000, 24, 2571−2588. (94) Alhumaizi, K.; Henda, R.; Soliman, M. Numerical analysis of a reaction-diffusion-convection system. Comput. Chem. Eng. 2003, 27, 579−594. (95) Nigam, K. M.; Nigam, K. D. P. Application of Galerkin technique to diffusion problems in tubular reactor. Chem. Eng. Sci. 1980, 35, 2358−2360. (96) Roussos, A. I.; Alexopoulos, A. H.; Kiparissides, C. Part III: Dynamic evolution of the particle size distribution in batch and continuous particulate processes: A Galerkin on finite elements approach. Chem. Eng. Sci. 2005, 60, 6998−7010. (97) Shafi, M. A.; Joshi, K.; Flumerfelt, R. W. Bubble size distribution in freely expanded polymer foams. Chem. Eng. Sci. 1997, 52, 635−644. (98) Solsvik, J.; Jakobsen, H. A. On the solution of the population balance equation for bubbly flows using the high-order least-squares method: Implementation issues. Rev Chem Eng 2013, 29, 63−98. (99) Coulaloglou, C. A.; Tavlarides, L. L. Description of interaction processes in agitated liquid−liquid dispersions. Chem. Eng. Sci. 1977, 32, 1289−1297. (100) Prince, M. J.; Blanch, H. W. Bubble coalescence and break-up in air-sparged bubble columns. AIChE J. 1990, 36, 1485−1499. (101) Morel, C.; Ruyer, P.; N, S.; Lavièville, J. M. Comparison of several models for multi-size bubble flows on an adiabatic experiment. Int J Multiphase Flow 2010, 36, 25−39. (102) Proot, M. M. J. The least-squares spectral element method. Ph.D. Thesis, Delt University of Technology, The Netherlands, 2003. (103) Jiang, B.-N. The least-squares Finite Element Method: Theory and Applications in Computational Fluid Dynamics and Electromagnetics; Springer: Berlin, 1998. (104) Bochev, P. B.; Gunzburger, M. D. Finite element methods of least-squares type. SIAM Rev. 1998, 40, 789−837. (105) Bochev, P. B.; Gunzburger, M. D. Least-squares finite-element methods for optimization and control problems for the Stokes equations. Comput. Math. Appl. 2004, 48, 1035−1057. 16012



Article

(130) Lasheras, J. C.; Eastwood, C.; Martínez-Bazán, C.; Montañeś , J. L. A review of statistical models for the break-up of an immiscible fluid immersed into fully developed turbulent flow. Int J Multiphase Flow 2002, 28, 247−278. (131) Ferziger, J. H.; Peric̀, M. Computational Methods for Fluid Dynamics, 3rd ed.; Berlin: Springer, 2002. (132) Patruno, L. E.; Marchetti, J. M.; Jakobsen, H. A.; Svendsen, H. F. Modeling of droplet entrainment in annular two-fluid three-phase dispersed flow. The 13th International Topical Meeting on Nuclear Reactor Thermal Hydraulics (NURETH-13), N13P1056, 2009. (133) Solsvik, J.; Tangen, S.; Jakobsen, H. A. Evaluation of weighted residual methods for the solution of the pellet equations: The orthogonal collocation, Galerkin, tau and least-squares methods. Comput. Chem. Eng. 2013, 58, 223−259. (134) Rout, K. R. A study of the sorption enhanced steam methane reforming process. Ph.D. Thesis, Norwegian University of Science and Technology (NTNU), 2012. (135) Szegö, G. Orthogonal Polynomials; American Mathematical Society: New York: , 1939. (136) Golub, G. H.; Welsch, J. H. Calculation of Gauss quadrature rules. Math. Comput. 1969, 23, 221−230. (137) Golub, G. H. Some modified matrix eigenvalue problems. SIAM Rev. 1973, 15, 318−334. (138) Crandall, S. H. Engineering Analysis; McGraw-Hill; New York, 1956. (139) Vichnevetsky, R. Use of functional approximation methods in the computer solution of initial value partial differential equation problems. IEEE Trans. Comput. C-18 1969, 18, 499−512.

(106) Bochev, P.; Gunzburger, M. D. Least-Squares Finite Element Methods; New York: Springer, 2009. (107) Bolton, P.; Thatcher, R. W. On mass conservation in leastsquares methods. J. Comput. Phys. 2005, 203, 287−304. (108) Pontaza, J. P. Least-squares variational principles and finite element method: theory, formulation, and models for solid and fluid mechanics. Ph.D. Thesis, Texas A&M University, 2003. (109) Pontaza, J. P. A least-squares finite element formulation for unstady incompressible flows with improved velocity-pressure coupling. J. Comput. Phys. 2006, 217, 563−588. (110) Pontaza, J. P. A new consistent splitting scheme for incrompressible Navier-Stokes flows: A least-squares spectral element implementation. J. Comput. Phys. 2006, 225, 1590−1602. (111) Pontaza, J. P.; Reddy, J. N. Spectral/hp least-squares finite element formulation for the incompressible Navier-Stokes equation. J. Comput. Phys. 2003, 190, 523−549. (112) Pontaza, J. P.; Reddy, J. N. Space-time coupled spectral/hp least squares finite element formulation for the incompressible NavierStokes equation. J. Comput. Phys. 2004, 190, 418−459. (113) Pontaza, J. P.; Reddy, J. N. Least-squares finite element formulations for viscous incompressible and compressible fluid flows. Appl. Mech. Eng. 2006, 195, 2454−2494. (114) Proot, M. M. J.; Gerritsma, M. I. A least-squres spectral element formulation for Stokes problem. J. Sci. Commut. 2002, 17, 285−296. (115) Proot, M. M. J.; Gerritsma, M. I. Mass- and momentum conservation of the least-squares spectral element method for the Stokes problem. J. Sci. Commut. 2006, 27, 389−401. (116) Solsvik, J.; Jakobsen, H. A. Effect of Jacobi polynomials on the numerical solution of the pellet equation using the orthogonal collocation, Galerkin, tau and least squares methods. Comput. Chem. Eng. 2012, 39, 1−21. (117) Ritz, W. Uber eine neue Methode zur Lsnung gewisses Variationsprobleme der mathematishen Physik. J. Reine. Angew. Math. 1908, 135, 1−61. (118) Rayleigh, J. W. In finding the correction for the open end of an organ-pipe. Phil. Trans. 1870, 161, 77. (119) Galerkin, B. G. Series solution of some problems in elastic equilibrium of rods and plates. Vestn. Inzh. Tech. 1915, 19, 897−908. (120) Villadsen, J. V.; Steward, W. E. Solution of boundary-value problems by orthogonal collocation. Chem. Eng. Sci. 1967, 22, 1482− 1501. (121) Villadsen, J. Selected approximation methods for chemical engineering problems; København: Danmarks Tekniske Højskole, 1970. (122) Villadsen, J.; Michelsen, M. L. Solution of Differential Equation Models by Polynomial Approximation; Prentice-Hall: Englewood Cliffs, NJ, 1978. (123) Michelsen, M. L.; Villadsen, J. Polynomial solution of differential equations. 1981; In Foundations of Computer-Aided Chemical Process Design; Mah, R S H and Seider, W D, Eds.; Engineering Foundation: New York, 341−368. (124) Finlayson, B. A. The Method of Weighted Residuals and Variational Principles; Mathematics in Science and Engineering; New York: Academic Press, 1972; Vol. 87. (125) Rice, R. G.; Do, D. D. Applied Mathematics and Modeling for Chemical Engineers; John Wiley & Sons: New York, 1995. (126) Grienberger, J. Untersuchung und Modellierund von Blasensäulen. Ph.D. Thesis, Der Technishen Fakultät der Universität Erlangen-Nürnberg, 1992. (127) Kostoglou, M.; Karabelas, A. J. On sectional techniques for the solution of the breakage equation. Comput. Chem. Eng. 2009, 33, 112− 121. (128) Liao, Y.; Lucas, D. A literature review of theoretical models for drop and bubble breakup in turbulent dispersions. Chem. Eng. Sci. 2009, 64, 3389−3406. (129) Liao, Y.; Lucas, D. A literature review on mechanisms and models for the coalescence process of fluid particles. Chem. Eng. Sci. 2010, 65, 2851−2864. 16013


Evaluation of Weighted Residual Methods for the ... - ACS Publications

Recommend Documents