Variational iteration method and projection method solution of the spatially distributed population balance equation

Abstract In this work, two major hydrodynamic parameters, the holdup of the dispersed phase and the Sauter diameter, are considered. This is done for describing the hydrodynamics of interacting liquid–liquid dispersions using different particle breakage, coalescence and growth models in a particle population balance model. Based on the semi-analytical solution method of the population balance, namely, the variational iteration method (VIM), different process cases have been performed, and it is possible to find the exact solution or a closed approximate solution of a problem. For the simultaneous growth and coalescence terms comparisons between the present method and projection method which include discontinuous Galerkin and collocation techniques are made, respectively. The VIM technique overcomes the difficulties of discretization of the variables, introduces an efficient algorithm that improves the standard discretization method and is able to handle quite successful these process of population balance equations. The results are encouraging and the new method has proven to be suitable to predict holdup and Sauter diameter profiles.


Introduction
The determination of the total number dispersed phase holdup and the particle size distributions in biphasic contactors is a key of chemical industry engineering processes such as drops or bubbles columns, dispersed phase polymerization and reactions in organic chemistry. Such systems have been modelled using population balance modelling, a method that describes the variations in size distribution of the dispersed phase as an averaged function of the behaviour of individual particles, drops or bubbles. This approach requires description of the interactions of the dispersed phase such as breakage, aggregation and growth.
The VIM was first proposed by He (1997), and it is well known that VIM provides the most versatile tools available for nonlinear analysis problems (He, 2000(He, , 2004(He, , 2006(He, , 2007Mohyud-Din, Sikander, & Naveed Ahmed, 2016). This method was shown to be effective, easy and accurate in solving a large class of nonlinear problems. Generally, one or two iterations lead to highly accurate solutions. This method is, in fact, a modification of the general Lagrange multiplier method into an iteration method, which is called correction functional.
The fixed-pivot technique of Kumar and Ramakrishna (1996) was used in the simulation of a bubble column through a mixed Euler-Lagrangian formulation (Campos & Lage, 2003) and the authors tried to solve it for droplet breakage and growth using the successive generation method (Liou et al., 1997). Unfortunately, they were not able to obtain an analytical solution in the closed from, and hence their semi-analytical solution requires further sophisticated programming, as they pointed out. Attarakih et al. (2004) developed a numerical algorithm to solve the PBE describing the hydrodynamics of interacting liquid-liquid dispersions showing both spatial and time dependencies as well as droplet interactions (breakage and coalescence), and this is accomplished by generalizing the fixed-pivot.
This work developed a semi-analytical methodology to solve the dynamic multidimensional PBE including advection in the internal and external coordinates, breakage, aggregation and growth, based on the VIM and also solve the PBE for the growth process by the discontinuous Galerkin method, which first was applied by Sandu (2004) for the aerosol dynamics. The discontinuous Galerkin method was designed as an effective numerical methods for solving hyperbolic conservation laws, which may have discontinuous solutions (Cockburn, 2003;Cockburn, Lin, & Shu, 1989).

Model equations
In the present model the basic variable is a drop size distribution function P(t,z,v) that represents the volume fraction of drops with volume v in an unit volume of the column at level z and time t (Attarakih et al., 2004). The local holdup of the dispersed phase can be calculated from P(t,z,v) as follows: The volume balance equation of unit volume of the column can be expressed as: (2) with a transient and a convective term, where v d is the dispersed phase velocity, balanced against a back mixing term expressed with the dispersion coefficient D ax . Q d and A are the dispersed phase flow and the column cross-sectional area, respectively. The solvent feed entering at the level z d of the column is handled as a point source by Dirac's d-function (Kronberger, Ortner, Zulehner, & Bart, 1995). The break-up, coalescence and growth process of drops are taken into account in the term p v , which is given by: The first integral accounts for gain and loss due to break-up of the mother particle u according to the daughter particle distribution b. The second integral holds for similarly for aggregation according to the aggregation frequency x.
In order to apply the proposed method we reformulate Eq. (2) as follows: where x ¼ z À z d and the boundary condition is given by Moreover and by assuming the equation of motion with uniform particle velocity; negligible diffusion flux, and making the necessary variable transformation using the chain rule, Eq. (4) could be reduced to (Attarakih et al., 2004): with h define the relative time as The following is the list of relevant combinations of processes for which the continuous PBE has been solved analytically. Case study I. Pure breakage: Case study II. Pure aggregation: Case study III. Pure growth: Case study IV. Aggregation and growth: Case study V. Breakage and growth: Case study VI. Simultaneous breakage, aggregation and growth:

Variational iteration method
As a first step to obtain the solutions for the above set of PBEs using the VIM, we consider the following functional equation: where L is a linear operator, N a nonlinear operator and g h; x; v ð Þ a source term. According to the VIM, a functional correction can be constructed as follows: (15) where k n ð Þ is a general Lagrangian multiplier which can be identified optimally via variational theory and e p n is a restrictive variation meaning f @p n , He (1997). Therefore, the solution is given by:

Projection method
The objective in this section is to develop the finite element expansion coefficients scheme based on discontinuous Galerkin and collocation methods for the solution of the simultaneous growth and coalescence process, respectively.

Discontinuous Galerkin method for the growth equation
To introduce the basic idea of the discontinuous Galerkin method, consider the following initial value problem: with the initial condition, Eq. (5). This has the form of conservation law with flux function f and source terms It is assumed that the spatial domain X is periodic and partitioned into nonoverlapping intervals I j; j ¼ 1, … , N j . The centre of bin I j is denoted by v j .
A weak formulation of the problem is obtained by multiplying (18) by an arbitrary smooth test function u v ð Þ and integrating over an interval I j , ð Integrating the second term of (18) by parts yields where are the values of the function u v ð Þ at the end points v jþ1 and v j of the element I j , respectively.
At an interface between elements (e.g. the points v jþ1 and v j Þ , the flux function f is not uniquely defined, and a suitable numerical flux must be determined according to the classical finite-volume method. For example, the nonlinear flux function Á that depends on two values, the left and right limits of the discontinuous function p evaluated at the interface v jþ1 such thatf ðf Þ jþ1 x ð Þ 1 . However, the Godunov numerical flux was chosen for the present study. For the approximate solution p h h; x; v ð Þ, the discontinuous Galerkin space discretization based on the weak formulation (19) is written as follows (Cockburn et al., 1989): We consider the approximations solutions p h for each interval I j ; p h ðxÞ j | is a polynomial of degree k. We take as the local basis function the suitably scaled Legendre polynomials, that is, for and P l is the lth Legendre polynomial. Since these polynomials are orthogonal, that is, , the mass matrix is diagonal. Indeed, the weak formulation (21) takes the following simple form: For each interval I j and each l ¼ 0, … , k, we have A high-order Gaussian quadrature rule is applied to evaluate the integral on the right-hand side of Eq. (22).
In order to obtain the approximation of the Eq. (6) which is solved in finite-volume interval I j by the collocation method for the break-up and coalescence, the approximate solution is then obtained by inserting the function p h h; x; v ð Þ into the weak formulation of (6) and multiplying both sides of the aggregation and break-up equation by choosing the test function as Dirac functions h j ðvÞ ¼ @ðv À v c j Þ with v c j the collocation points. In this section, we briefly display the application of these discretizations. For more details the reader should consult references Sandu, 2004;Sandu & Borden, 2003).

Application and results
Let us consider the spatially distributed one-dimensional PBE Eq. (6) subject to the boundary condition pðh; 0; vÞ ¼ ve Àv for the following cases: In the following, we solve this problem by the VIM.
Integrating Eq. (8) with respect to x we have Hence we calculate the general term as Then: which converges to the exact solution. Now rewriting this equation in terms of the original variables we get (29) Figure 1 shows the prediction of the number density by the VIM at z d ¼ 0.3 (inlet) and z ¼ 1.8 (outlet). It must be noted from these profiles that the outlet number density with linear breakage rate gives much production of particles compared with the inlet density as expected. column increases and the Sauter diameter decreases, as shown in Figures (2a) and (2b), respectively. As can be seen from the figure, the semi-analytical and numerical results of VIM and collocation, respectively, are very distinguishable. The solid line in Figure 3 shows the holdup versus column height obtained by the VIM method and the markers by the collocation approach.
Case 2. Aggregation with x v; u ð Þ ¼ constant and v d ¼ constant: À 3 4 e Àv v 2 x 2 þ 1 8 e Àv v 3 x 2 þ :::; (33) then we calculate the general term as The closed-form solution can be written as The solution in terms of the original variables is given by (36) Figure 4 shows the prediction of the number density by the VIM at z d = 0.3 (inlet) and z ¼ 2 (outlet). It must be noted from these profiles that the outlet number density with constant aggregation rate gives low production of drops compared with the inlet density, as expected. Figure 5 concerns the total number of particles; Figure   lower total number of particles and a higher Sauter diameter, as shown in Figures (5a) and (5b), respectively. Again, it can be seen from the figure that the semi-analytical and numerical results of VIM and collocation, respectively, are in good agreement. The solid line in Figure 6 shows the holdup versus column height obtained by the VIM method and the markers by the collocation approach.
Case 3. Growth with G v ð Þ ¼ v and v d ¼ constant: þ v e Àv À e Àv x þ 1 2 e Àv x 2 :::; then the general can be obtained as The solution can be written as which can be written as Figure 7 presents a comparison between the discontinuous Galerkin method and the VIM results for particle growth Eq. (10).
Case 5. Simultaneous breakage and growth with G v ð Þ ¼ v and v d ¼ constant: e Àv x 2 :::; then the general term is obtained as With the closed solution written as Then the solution is (63) Figure 9 shows a comparison between the discontinuous Galerkin method and VIM results for particle growth and breakage Eq. (12).
Case 6. Simultaneous breakage and coalescence and growth with G v ð Þ ¼ v , x v; u ð Þ ¼ constant and v d ¼ constant: þ 1 8 e Àv v 3 x 2 ::: (67) Figure 10 shows different cases for the total number of particles such as pure breakage, pure aggregation, pure growth, aggregation and growth, and simultaneous breakage and aggregation and growth are considered. The results can be regarded as a sensitivity analysis of the base concept of a population balance.

Constant aggregation rate
Linear aggregation rate Figure 11. Comparison between the total number of particles versus column height for the case of pure aggregation process with constant aggregation rate and linear aggregation rate.
þ v e Àv À e Àv x þ 3 4 e Àv x 2 À 1 6 e Àv x 3 þ v 5 1 480 e Àv x 2 À 1 720 e Àv x 3 þ v 3 e Àv x 12 À 1 8 e Àv x 2 þ 1 24 e Àv x 3 ; Using the above terms, the general term is deduced as Gamma 2n ð Þ e Àv ; The above series can be generalized as follows:  Figures 12 and 13 show the prediction of the volume density by the VIM at z d ¼ 0.3 (inlet) and z ¼ 2.2 (outlet) for the breakage and aggregation respectively with linear velocity. These profiles show that the outlet number densities with linear velocity give low production of drops for the aggregation and greater production of drops for the breakage compared with the inlet densities, as expected.

Conclusions
This paper explored the use of the VIM applied to hydrodynamics simulation based on particles population balance model for bubbles or droplets column. The comparison between the present method and the projection method which includes collocation (aggregation, breakage) and discontinuous Galerkin (growth) techniques is made. The results showed that the VIM eliminates complex calculations and provides highly accurate numerical solutions without spatial discretizations for the PBEs. It is also worth noting that the VIM methodology displays a fast convergence of the solutions. The illustrations show the rapid convergence of the solutions just as in a closed-form solutions.

Volume density
Outlet density at z 2.2 Inlet density at zd 0.3 Figure 13 Volume density for the pure aggregation with linear velocity by the variational iteration method at z d ¼ 0.3 (inlet) and z ¼ 2.2 (outlet).