A General Procedure for Solving the Population Balance Equation in Flocculating Suspensions

A General Procedure for Solving the Population Balance Equation in Flocculating Suspensions

Leong Y Yeow

The University of Melbourne

Jong-Leng Liow

University of NSW Canberra @ADFA

Yee-Kwong Leong

The University of Western Australia

In a flocculating suspension large particles are continuously being formed by agglomeration of smaller ones and at the same time small particles are being generated by fragmentation of larger ones.  At any time t the particle size distribution (PSD) is described by the number density function  which gives the number concentration density of size v particles. 

The evolution of the PSD n(t,v) is governed by the following population balance equation (PBE)

Eq (1)

In this equation, k_A(v_s,v_t) describes the kinetics of formation large particles [v_s+v_t→v_u] by aggregation, k_F(v_r,v_p) describes the kinetics of the production of small particles [v_r →v_p+v_q] by fragmentation.

The first integral on the RHS of Eq. (1) gives the rate of generation of size v particles as a result of the agglomeration of particles of size v_s with size v – v_s. The second integral is the rate of disappearance of size v particles as a result of the agglomeration with particles of size v_t. The third integral is the rate of generation of size v particles arising from the fragmentation of all size v_r particles into a size v particle and a second size v_r – v particle.  The final integral accounts for the rate of disappearance of size v particles as a result of their fragmentation into all possible combinations of two smaller particles.  Together these four terms give the net rate of change of the number concentration density n of size v particles.

The PBE is an integro-differential equation that does not have a simple analytical solution and for practical flocculation problems, solving it numerically results in a major difficulty because of the large span in particle volume; usually spread over 9 or more decades of particle volume or 3 decades of particle diameter.  Uniform discretization of the independent variable v then leads to inadequate numerical representation for small particle sizes.

A popular procedure for overcoming the numerical problem is to discretize the span of particle size geometrically so that Δv_j+1 = 2Δv_j. This rectifies the under representation of small particles but is less intuitive. In the present method this problem is surmounted by the introduction of the logarithmic variable s = log₁₀v as the new independent variable providing both flexibility as well as a continuous function for the particle size which can be easily be handled within the PBE. The new dependent variable is n(t,s) and the PBE is then reformulated. In the computation that follows the span of s is discretised into a uniform grid i.e. Δs_j+1 = Δs_j..  This ensures that the PBE is adequately represented over the entire span of particle volume.

Another key development is to apply trapezoidal rule to approximate the four integrals on the RHS of the PBE.  After discretization, the integro-differential equation becomes a set of first order ODE for the discretised analogue of n(t, s) i.e. n₁(t), n₂(t),… n_j(t)… n_N(t).  where N is the number of discretization steps.  This set of ODE can be written in a matrix form as

The matrix a depends on the aggregation and fragmentation kinetics.  For a given kinetics, a_ij are known functions of the discretised unknowns n_i(t).  This is a set of nonlinear first order ODEs in time that we solved numerically with n_i(0), the PSD at the start of the flocculation, as the initial condition. This method of handling the integro-differential is an extension of the methods of line (MoL) technique used to solve partial differential equations. 

To cope with wide span in particle volume, the number of discretization interval (N) can be as larger than 801. Although a is large matrix, it is sparse and can be handled efficiently by numerical procedures found in most scientific computing software.

 

Results

Typical computed results for different simulated kinetics and initial PSD are presented in Figures 1 to 3 in dimensionless form, where τ, V and N are the dimensionless time, particle volume and number concentration density function respectively. 

In addition, the results are also presented as volume fraction density function G(X) as a function of X= log₁₀V.  Experimental PSD are usually reported as G(X).  The area under each G(X) curve is by definition equal to unity and provides a means to validate the results. The area under the G(x) in all the cases remained below 2% confirming the reliability of the computed results.

For all the aggregation kinetic plots each curve represents a constant value of V_t, and all fragmentation kinetic plot each curve represent a constant value of V_p. In the number concentration density plots and the particulate volume plots the discrete points represent the initial PSD.  The curves show the evolving PSD with a constant dimensionless time on each curve. The results in Figure 3 for rain drops show that the method is capable of capturing the evolution of the PSD for a natural event, particularly for the small values of particles sizes.

    (a)                                (b)                            (c)                                  (d)

Figure 1 (a) and (b) Aggregation and fragmentation kinetics; (c) Dimensionless number concentration density function; (d) Volume fraction density function.

 

    (a)                                    (b)                                (c)                                      (d)

Figure 2 (a) and (b) Aggregation and fragmentation kinetics; (c) Dimensionless number concentration density function;(d) Volume fraction density function.

(a)                                              (b)                                          (c)

Figure 3 Coalescence of rain droplets under constant aggregation kinetics.  (a) Dimensionless number concentration density function; (b) Replot of (a) with log scale on both axes.  (c) Volume fraction density function.  Dark continuous curve are the computed results of MoL and the lighter curves are exact results of Scott (1968).

Reference

Scott WT. Analytic studies of cloud droplet coalescence I.  J. Atmos Sci. 1968; 25:54-65.