(199c) 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.