# Comparison of the Accuracy and Performance of Different Numbers of Classes in Discretised Solution Method for Population Balance Model

Authors: Zhenliang Li; Zhien Zhou; Sheng Zhang; Hongqiang Jiang
Journal: International Journal of Chemical Engineering
 (2016)

---

## Abstract

One way of solving population balance model (PBM) in a time efficient way is by means of discretisation of the population property of interest. A computational grid, for example,v
i
+
1
=
k
v
i (v
i is the volume of particle in class i), could be used to classify the particles in discretisation techniques. However, there are still disagreements in the appropriate number of classes divided by the grids. In this study, the different numbers of classes for solving PBM were compared in terms of accuracy and performance to describe the particle size distribution (PSD) from the flocculation of activated sludge. It is found that the simulated PSDs are similar to the experimental data for all the geometric grids (
v
i
+
1
:
v
i
≤
2
), and there is no obvious difference among the values of calibrated parameter, ratio of breakage rate coefficient and collision efficiency, for each velocity gradient. However, the simulation results with less error could be obtained with larger number of classes, and more computational times, which show exponential relationship with the number of classes, are needed. Considering numerical accuracy and efficiency, the classes 35 or a geometric grid with factor 1.6, aligning with the Fibonacci sequence (
v
i
+
v
i
-
1
≈
v
i
+
1
), is recommended for the particles in the size range of 5.5~1086 μm.


---

## Body

## 1. Introduction

Population balance model (PBM) has been widely used to describe the processes involving dynamical behaviour of population properties. Applications can be found in a variety of scientific areas such as flocculation, crystallisation, polymerisation, granulation, and aerosol and cell culture dynamics [1]. The numerical techniques to solve this type of equations are summarised by Ramkrishna [2]. One of these techniques resulting in an acceptable computational time and accuracy is the discretisation of the particle size.Discretisation techniques divide the property range of interest into a finite number of classes (n) according to a given computational grid, transforming the PBM into a set of n ordinary differential equations that can be solved simultaneously [3]. A computational grid, for example, v
i
+
1
=
k
v
i (v
i is the volume of particle in class i), could be used to classify the particles for solving PBMs. However, there are still disagreements in the appropriate number of classes. McAnally and Mehta [4] concluded that an available geometric grid used to classify the particles should be v
i
+
1
:
v
i
≤
2 in order to maintain the continuity of particles distribution after aggregation. Calabrese et al. [5] demonstrated that the geometric grid with factor 2 is not as suitable as another grid which kept in accordance with the Fibonacci sequence, v
i
+
1
=
v
i
+
v
i
-
1. Xu et al. [6] suggested a computational grid v
i
+
1
:
v
i
=
2
1
/
4 to achieve a reasonable numerical accuracy. Verney et al. [7] believed that the number of classes n
≥
14 could obtain the satisfying simulation results when the particle size is in the range of 4~800 μm, whereas Mietta et al. [8] recommended the number of classes n
≥
30 for nearly the same range of particles. Absolutely, the computational cost increases with the increasing number of classes. Acceptable computational speeds are important when using the model for parameter estimation or in combination with CFD [9]. Therefore, the problem needs to be firstly solved in the application of PBM: what is the appropriate number of classes in discretisation techniques? In addition, it is still not clear whether the number of classes has influence on the results of parameters calibration.In this study, different numbers of classes in the discretised solution for PBM were compared, in terms of accuracy and performance to describe the particle size distribution (PSD) from the flocculation of activated sludge, and the appropriate number of classes was recommended according to a reasonable compromise between accuracy and efficiency of PBM.


## 2. Theory and Experiments

### 2.1. Theoretical Section

The discretised PBM developed by Kumar and Ramkrishna [10] is used to describe the variation rate in the number of particles with a given size:(1)
d
N
i
d
t
=
∑
j
,
k
v
i
-
1
≤
v
j
+
v
k
≤
v
i
+
1
j
≥
k
1
-
1
2
δ
j
,
k
η
i
α
β
j
,
k
N
j
N
k
-
N
i
∑
k
α
β
i
,
k
N
k
+
∑
j
≥
i
γ
j
,
i
S
j
N
j
-
S
i
N
i
,where N
i is the number concentration of particles of class i, δ
j
,
k is Dirac delta function, α is the collision efficiency, β
(
j
,
k
) is the collision frequency of particles in classes j and k, S
(
i
) is the breakage rate of particles in class i, γ
j
,
k is the breakage distribution function defining the fraction of daughter-particle of class size i breaking from a particle of class size j, and η
i is a proportional coefficient assigning the fraction of the particle v
i from the aggregate (v
j
+
v
k):(2)
η
i
=
v
i
+
1
-
v
j
+
v
k
v
i
+
1
-
v
i
v
i
≤
v
j
+
v
k
≤
v
i
+
1
v
j
+
v
k
-
v
i
-
1
v
i
-
v
i
-
1
v
i
-
1
≤
v
j
+
v
k
≤
v
i
.The collision efficiency is generally considered as a constant which needs to be calibrated with the experiments [11]. The collision frequency due to shear rate has usually been considered as function of particles size L, as well as the breakage rate [12, 13]:(3)
β
i
,
j
=
G
6
L
i
+
L
j
3
,where G is the average velocity gradient:(4)
G
=
ε
ν
1
/
2
,where ε represents the homogeneous turbulent energy dissipation rate of the mixed tank and ν is the kinematic viscosity of the suspending fluid.(5)
S
i
=
E
L
i
,where E is the breakage rate coefficient.For a geometric grid adopted to solve the population balance model, the binomial breakage function can be used for the daughter-particle distribution [14]:(6)
γ
j
,
i
=
j
i
p
i
1
-
p
j
-
i
,
i
=
1,2
,
…
,
j
,
p
=
i
-
C
p
i
,where C
p is the parameter which is defined to determine the location at which the probability mass function reaches its peak and can be calculated from(7)
C
p
=
log
⁡
2
log
⁡
k
.


### 2.2. Modelling Method

Integrating (1) involved computing the discrete Δ
N
i
/
Δ
t of d
N
i
/
d
t values. The geometric grids with a variable factor of k were used to classify the particle according to the size. The solution equations were derived using the Euler method, which involved establishing 0.1~1 s as the iterative calculation of Δ
t to maintain the calculation stability [15]. The final state of flocculation is that particle aggregation and breakage counterbalance each other; that is, d
N
/
d
t = 0 in (1). Connecting (5) and (1) yields(8)
d
N
i
d
t
=
α
∑
j
,
k
v
i
-
1
≤
v
j
+
v
k
≤
v
i
+
1
j
≥
k
1
-
1
2
δ
j
,
k
η
i
β
j
,
k
N
j
N
k
-
N
i
∑
k
β
i
,
k
N
k
-
E
L
i
N
i
-
∑
j
≥
i
γ
j
,
i
L
j
N
j
=
0
.Equation (8) shows the correlation between the breakage rate coefficient E and the collision efficiency α. Given the ratio of breakage rate coefficient and collision efficiency E
/
α, the number concentration of particles in each class could be obtained. The volume percentage distribution was chosen as the fitting variables of the minimum error for parameters calibration [16]:(9)
E
r
r
=
1
n
∑
i
=
1
n
m
i
-
m
′
i
,where m
(
i
) and m
′
(
i
) were the measured and simulated volume percentage of particles in class i, respectively.The optimization results of parameter, the ratio of breakage rate coefficient and collision efficiency, were estimated by the exhaust algorithm in the range of possible values. The computational speed was recorded during the simulation using Matlab 2012 (MathWorks Inc., USA) on the PC with Intel i5 processor systems (3.0 GHz).


### 2.3. Experimental Section

Flocculation experiments were conducted in a mixing tank (ZR4-6, China) using activated sludge collected from a municipal wastewater treatment plant. The volume of mixed liquid in the tank was 1 L and the concentration of sludge was 0.1 kg·m−3. The PSD (or volume percentage distribution) of activated sludge during flocculation was measured using a laser particle size analyzer (S3500, Microtrac, in the range of 0.02~2800 μm) at the velocity gradients of 28.2 s−1, 64.7 s−1, and 101.7 s−1, which are the typical values of velocity gradient for flocculation occurrence. In order to reduce the possible effect on the particle size distribution caused by the difference of flowing shear between the sampling pipe and the mixing tank, the velocity gradients in the sampling pipe were set to be similar to those in the mixing tank by controlling suitable sample flow rates in the sampling pipe [17].The raw output of the laser particle size analyzer is volume percentage distribution for a given grid. However, PBM is formulated on a number concentration basis and might use a different grid. The cumulative volume percentage distribution was first calculated from the raw volume percentage distribution and then interpolated at the pivots of the new grid, allowing the recalculation of the volume percentage distribution [18]. Finally, the initial number concentrations of the models were generated by transforming the volume distribution data to a number distribution.


## 2.1. Theoretical Section

The discretised PBM developed by Kumar and Ramkrishna [10] is used to describe the variation rate in the number of particles with a given size:(1)
d
N
i
d
t
=
∑
j
,
k
v
i
-
1
≤
v
j
+
v
k
≤
v
i
+
1
j
≥
k
1
-
1
2
δ
j
,
k
η
i
α
β
j
,
k
N
j
N
k
-
N
i
∑
k
α
β
i
,
k
N
k
+
∑
j
≥
i
γ
j
,
i
S
j
N
j
-
S
i
N
i
,where N
i is the number concentration of particles of class i, δ
j
,
k is Dirac delta function, α is the collision efficiency, β
(
j
,
k
) is the collision frequency of particles in classes j and k, S
(
i
) is the breakage rate of particles in class i, γ
j
,
k is the breakage distribution function defining the fraction of daughter-particle of class size i breaking from a particle of class size j, and η
i is a proportional coefficient assigning the fraction of the particle v
i from the aggregate (v
j
+
v
k):(2)
η
i
=
v
i
+
1
-
v
j
+
v
k
v
i
+
1
-
v
i
v
i
≤
v
j
+
v
k
≤
v
i
+
1
v
j
+
v
k
-
v
i
-
1
v
i
-
v
i
-
1
v
i
-
1
≤
v
j
+
v
k
≤
v
i
.The collision efficiency is generally considered as a constant which needs to be calibrated with the experiments [11]. The collision frequency due to shear rate has usually been considered as function of particles size L, as well as the breakage rate [12, 13]:(3)
β
i
,
j
=
G
6
L
i
+
L
j
3
,where G is the average velocity gradient:(4)
G
=
ε
ν
1
/
2
,where ε represents the homogeneous turbulent energy dissipation rate of the mixed tank and ν is the kinematic viscosity of the suspending fluid.(5)
S
i
=
E
L
i
,where E is the breakage rate coefficient.For a geometric grid adopted to solve the population balance model, the binomial breakage function can be used for the daughter-particle distribution [14]:(6)
γ
j
,
i
=
j
i
p
i
1
-
p
j
-
i
,
i
=
1,2
,
…
,
j
,
p
=
i
-
C
p
i
,where C
p is the parameter which is defined to determine the location at which the probability mass function reaches its peak and can be calculated from(7)
C
p
=
log
⁡
2
log
⁡
k
.


## 2.2. Modelling Method

Integrating (1) involved computing the discrete Δ
N
i
/
Δ
t of d
N
i
/
d
t values. The geometric grids with a variable factor of k were used to classify the particle according to the size. The solution equations were derived using the Euler method, which involved establishing 0.1~1 s as the iterative calculation of Δ
t to maintain the calculation stability [15]. The final state of flocculation is that particle aggregation and breakage counterbalance each other; that is, d
N
/
d
t = 0 in (1). Connecting (5) and (1) yields(8)
d
N
i
d
t
=
α
∑
j
,
k
v
i
-
1
≤
v
j
+
v
k
≤
v
i
+
1
j
≥
k
1
-
1
2
δ
j
,
k
η
i
β
j
,
k
N
j
N
k
-
N
i
∑
k
β
i
,
k
N
k
-
E
L
i
N
i
-
∑
j
≥
i
γ
j
,
i
L
j
N
j
=
0
.Equation (8) shows the correlation between the breakage rate coefficient E and the collision efficiency α. Given the ratio of breakage rate coefficient and collision efficiency E
/
α, the number concentration of particles in each class could be obtained. The volume percentage distribution was chosen as the fitting variables of the minimum error for parameters calibration [16]:(9)
E
r
r
=
1
n
∑
i
=
1
n
m
i
-
m
′
i
,where m
(
i
) and m
′
(
i
) were the measured and simulated volume percentage of particles in class i, respectively.The optimization results of parameter, the ratio of breakage rate coefficient and collision efficiency, were estimated by the exhaust algorithm in the range of possible values. The computational speed was recorded during the simulation using Matlab 2012 (MathWorks Inc., USA) on the PC with Intel i5 processor systems (3.0 GHz).


## 2.3. Experimental Section

Flocculation experiments were conducted in a mixing tank (ZR4-6, China) using activated sludge collected from a municipal wastewater treatment plant. The volume of mixed liquid in the tank was 1 L and the concentration of sludge was 0.1 kg·m−3. The PSD (or volume percentage distribution) of activated sludge during flocculation was measured using a laser particle size analyzer (S3500, Microtrac, in the range of 0.02~2800 μm) at the velocity gradients of 28.2 s−1, 64.7 s−1, and 101.7 s−1, which are the typical values of velocity gradient for flocculation occurrence. In order to reduce the possible effect on the particle size distribution caused by the difference of flowing shear between the sampling pipe and the mixing tank, the velocity gradients in the sampling pipe were set to be similar to those in the mixing tank by controlling suitable sample flow rates in the sampling pipe [17].The raw output of the laser particle size analyzer is volume percentage distribution for a given grid. However, PBM is formulated on a number concentration basis and might use a different grid. The cumulative volume percentage distribution was first calculated from the raw volume percentage distribution and then interpolated at the pivots of the new grid, allowing the recalculation of the volume percentage distribution [18]. Finally, the initial number concentrations of the models were generated by transforming the volume distribution data to a number distribution.


## 3. Results and Discussion

Figure1 shows the variation rate in the number of particles for three numbers of classes from different computational grids (take, e.g., the simulation of flocculation at the velocity gradient of 28.2 s−1). For t
=
300 s, d
N
/
d
t of some particles, especially whose sizes are smaller than 10 μm, is so far from zero (-
1.0
×
10
7~−2.0 × 106). For t
=
1800 s, d
N
/
d
t of all particles is nearly zero, indicating that the PSD at steady state of flocculation would be achieved soon. When t
≥
5400 s, d
N
/
d
t is far smaller than zero and unchanged with time; so 5400 s can be used for the total simulation time at which the steady-state PSD is obtained. Moreover, it seems that three different geometric grids make negligible difference to the variation rate in number of particles.Figure 1
Variation rate in number of particles by using three different computational grids for simulation of flocculation at velocity gradient of 28.2 s−1 (the ratio of breakage rate coefficient and collision efficiency E
/
α for each computational grid was listed in Table 1, and collision efficiency was 0.5 for all grids).
(a)
t
=
100 s
(b)
t
=
300 s
(c)
t
=
900 s
(d)
t
=
1800 s
(e)
t
=
5400 s
(f)
t
=
9000 sTable 1
Calibration parameters and error of simulation using different numbers of classes in discretised solution method for PBM at different velocity gradients.
Number of classes
Ratio of breakage rate coefficient and collision efficiency,E
/
α (m−1⋅s−1)
Average error of volumep
e
r
c
e
n
t
a
g
e
∗
Computational time (s)
Dimensionless number∗
∗, D
Note
G
=
28.2 s−1
24 (k= 2.0)
1032
0.0247
13.8
0.0900
The total simulation time at which the steady-state PSD can be obtained was 5400 s. The possible values of the parameterE
/
α were chosen in the range of [
900
,
1200
], and a constant step-length 1 was applied in the exhaust algorithm for parameter calibration.
28 (k= 1.8)
1045
0.0235
18.7
0.0980
35 (k= 1.6)
1025
0.0181
30.4
0.0760
48 (k= 1.4)
1052
0.0163
61.5
0.0977
88 (k= 1.2)
1050
0.0128
213.3
0.0900
G
=
64.7 s−1
24 (k= 2.0)
2602
0.0243
7.6
0.0900
The total simulation time at which the steady-state PSD can be obtained was 3000 s. The possible values of the parameterE
/
α were chosen in the range of [
2500
,
2700
], and a constant step-length 1 was applied in the exhaust algorithm for parameter calibration.
28 (k= 1.8)
2598
0.0231
10.1
0.1020
35 (k= 1.6)
2617
0.0165
14.6
0.0640
48 (k= 1.4)
2611
0.0150
25.5
0.0770
88 (k= 1.2)
2626
0.0123
89.4
0.0900
G
=
101.7 s−1
24 (k= 2.0)
4267
0.0287
5.5
0.0900
The total simulation time at which the steady-state PSD can be obtained was 2500 s. The possible values of the parameterE
/
α were chosen in the range of [
4100
,
4200
], and a constant step-length 1 was applied in the exhaust algorithm for parameter calibration.
28 (k= 1.8)
4185
0.0268
8.2
0.1130
35 (k= 1.6)
4228
0.0214
10.3
0.0842
48 (k= 1.4)
4226
0.0185
18.7
0.0959
88 (k= 1.2)
4199
0.0149
55.6
0.0900
∗: the cumulative volume percentage distribution was first calculated from the simulated PSD using different numbers of classes and then interpolated at the pivots of the measured PSD grid. Finally, the average error of volume percentage was recalculated by comparing with the measured PSD.
∗
∗: a dimensionless number; D: relating error and computational time; see (10).Similarly, it is concluded that 3000 s or 2500 s can be used as the total simulation time for obtaining the steady-state PSDs when the velocity gradient is 64.7 s−1 or 101.7 s−1, respectively.The simulation results of steady-state PSD are similar to experimental data for all the computation grids (Figure2), demonstrating that the geometric grids, v
i
+
1
:
v
i
≤
2, are available in discretised solution method for PBM. More number of classes significantly improves the simulation result of volume percentage distribution of small particles (size less than about 100 μm). However, the results are not satisfactory for large particles (size more than about 260 μm), and the deviation between the simulation results and experimental data increases with the increasing number of classes. This deviation might be due to the fact that the breakage rate of large flocs is overestimated according to (5), in which the breakage rate is simply proportional to floc size.Figure 2
Comparison of simulation results of volume percentage distribution by using different numbers of classes in discretised solution method for PBM at velocity gradient of 28.2 s−1 (the simulation results were recalculated from the cumulative volume percentage, which was first calculated from the volume percentage distribution simulated for different computation grids and then interpolated at the pivots of the measured grid).Table1 shows the average error of volume percentage between the simulation results and experimental data, which decreases with the increasing number of classes. However, more computational time for larger number of classes is needed. Comparing with the linear function, the exponential function gives better fitting results for the relationship between the average error of volume percentage and the number of classes (Figure 3) and for the relationship between the computational time and the number of classes (Figure 4).Figure 3
Relationship between average error of volume percentage and number of classes in discretised solution method for PBM.Figure 4
Relationship between computational time and number of classes in discretised solution method for PBM.A dimensionless number,D, was introduced to describe the coupled effects of error and computational time:(10)
D
=
E
×
T
,where E is the average error of volume percentage normalized in the range of [
0.1
,
0.9
] and T is the computational time normalized in the range of [
0.1
,
0.9
].According to the minimum value ofD (Table 1), the classes 35 or a geometric grid with factor 1.6, aligning with the Fibonacci sequence (v
i
+
v
i
-
1
≈
v
i
+
1), is recommended for the particles in the size range of 5.5~1086 μm.In Table1, the estimated values of parameter E
/
α are similar to each other for the same velocity gradient. Spearman’s rank correlation coefficient r
s was used to identify the difference among the estimated values of parameter E
/
α for the different numbers of classes:(11)
r
s
=
1
-
6
∑
d
i
2
m
m
2
-
1
,where d
i is the difference between the two ranks of each observation and m is the number of observations.r
s were calculated to be 0.6, 0.8, and 0.4 for the velocity gradients of 28.2 s−1, 64.7 s−1, and 101.7 s−1, respectively. They all are less than r
s
0.05
=
0.9. So it can be concluded that there is no obvious difference among the values of calibrated parameter using the different numbers of classes for solving PBM.


## 4. Conclusions

The different numbers of classes divided by different geometric grids for solving the PBM were compared in terms of accuracy and performance to describe the PSD from the flocculation of activated sludge. For a given simulated time (e.g., 5400 s), the steady-state PSDs could be obtained for three different geometric grids, and these grids make negligible difference to the variation rate in number of particles. The simulated PSDs are similar to the experimental data for the geometric grids (v
i
+
1
:
v
i
≤
2). However, more accuracy results could be obtained with larger number of classes, and higher computational cost is needed. The results indicate that the computational time is exponentially related to the number of classes. Considering numerical accuracy and efficiency, the classes 35 or a geometric grid with factor 1.6, aligning with the Fibonacci sequence (v
i
+
v
i
-
1
≈
v
i
+
1), is recommended for the particles in the size range of 5.5~1086 μm. In addition, it can be concluded that there is no obvious difference among the values of calibrated parameter, ratio of breakage rate coefficient and collision efficiency, for the simulations using the different numbers of classes. These results will be helpful to the application of PBM.



---

*Source: 8743105-2016-11-22.xml*