Development of a Population Balance Model of a pharmaceutical drying process and testing of solution methods

Computers and Chemical Engineering 50 (2013) 39–53

Contents lists available at SciVerse ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

Development of a Population Balance Model of a pharmaceutical drying process and testing of solution methods Séverine Thérèse F.C. Mortier a,b,2 , Krist V. Gernaey c , Thomas De Beer b,1 , Ingmar Nopens a,∗ a

BIOMATH, Department of Mathematical Modelling, Statistics and Bioinformatics, Faculty of Bioscience Engineering, Ghent University, Coupure Links 653, 9000 Ghent, Belgium Laboratory of Pharmaceutical Process Analytical Technology, Department of Pharmaceutical Analysis, Faculty of Pharmaceutical Sciences, Ghent University, Harelbekestraat 72, 9000 Ghent, Belgium c Center for Process Engineering and Technology, Department of Chemical and Biochemical Engineering, Technical University of Denmark, Building 229, 2800 Kgs. Lyngby, Denmark b

a r t i c l e

i n f o

Article history: Received 31 August 2012 Received in revised form 15 November 2012 Accepted 17 November 2012 Available online 23 November 2012 Keywords: Drying Pharmaceuticals Mathematical modelling PBM Solution methodology

a b s t r a c t Drying is frequently used in the production of pharmaceutical tablets. Simulation-based control strategy development for such a drying process requires a detailed model. First, the drying of wet granules is modelled using a Population Balance Model. A growth term based on a reduced model was used, which describes the decrease of the moisture content, to follow the moisture content distribution for a batch of granules. Secondly, different solution methods for solving the PBM are compared. The effect of grid size (discretization methods) is analyzed in terms of accuracy and calculation time. All tested methods are compared based on their ability to predict moment dynamics and the distribution, and their computational burden. The Method of Characteristics, a fast method, is able to calculate the distribution accurately with a coarse grid. The Quadrature Method of Moments requires even less calculation time, but results in a set of moments. © 2012 Elsevier Ltd. All rights reserved.

1. Introduction A fluidized bed drying process is a commonly used technique for the drying of pharmaceutical granules. Such a drying step is performed (Mortier et al., 2011) when a preceding wet granulation technique is used during the production of tablets. The drying step might influence the properties of the granules, so the choice of the drying technique is important (Hegedus & Pintye-Hódi, 2007). Following the drying step, the dry granules may be further processed in a milling step and/or a post-blending step before reaching the final tabletting step (Mortier, De Beer, et al., 2012). To ensure the quality of the tablets, the complete process should be controlled tightly. Nowadays, the emergence of continuous production processes enhances the necessity of developing mechanistic models. These models are in general advantageous in terms of process understanding and process control. Indeed, when developing a mechanistic model, the process knowledge about the most

∗ Corresponding author. E-mail addresses: [email protected] (S.T.F.C. Mortier), [email protected] (K.V. Gernaey), [email protected] (T. De Beer), [email protected] (I. Nopens). URL: http://biomath.ugent.be (S.T.F.C. Mortier). 1 Shared last authorship. 2 Tel.: +32 (0)9 264 61 96; fax: +32 (0)9 264 62 20. 0098-1354/$ – see front matter © 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compchemeng.2012.11.005

important input–output relations is coded in equations, and can then later on be helpful in detecting input variables which have a significant influence on the process. Most important, indeed, is that such a mechanistic model is actively used for in silico testing of a set of potential process operating strategies, e.g. by comparing different control strategies in a series of dynamic simulations, without the need of disturbing the actual full-scale process operation (Gernaey, Cervera-Padrell, & Woodley, 2012). When using this information to establish control loops with sufficient control authority, the continuous process can subsequently be controlled based on in-process measurements and real-time adjustment of input variables in order to stay inside the Design Space (Mortier et al., 2011). In previous work, a mechanistic model to describe the drying behaviour of single pharmaceutical granules was developed (Mortier, De Beer, et al., 2012). The wet pharmaceutical granules consist of a wet porous particle surrounded by a water layer at the surface. Hence, the model consists of two drying phases; water from the droplet evaporates in the first drying phase, characterized by a high evaporation rate. When the radius of the wet granule equals the radius of the dry granule, the second drying phase initiates, in which water evaporates at the interface between the wet core and the dry crust. The produced vapor diffuses through the pores until it reaches the ambient air, where it is removed through advection by the air flow. In the second phase the evaporation rate is much

40

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53

0.08 0.07 0.06

X (%)

0.05 0.04 0.03 0.02 0.01 0

0

100

200

300 t (s)

400

500

600

Fig. 1. Evolution of the moisture content for one single pharmaceutical granule (Tg = 55 ◦ C, Vg = 200 m3 /h, Rp = 6e − 4 m.

lower (Mortier, De Beer, et al., 2012). The evolution of the moisture content for one single pharmaceutical granule is presented in Fig. 1. This drying model was calibrated and validated using experimental data (Fig. 2). Subsequently, a model reduction was conducted (Fig. 2). It consists of a Global Sensitivity Analysis (GSA) to detect the most sensitive degrees of freedom followed by a model reduction. The most important degrees of freedom are the temperature (Tg ) and velocity (Vg ) of the fluidized bed incoming gas and the dry particle radius (Rp ). The outcome of the model reduction was an empirical model (considerably less complex than the original model) able to describe the decrease of the wet radius (Rw ), which is a representation of the moisture content, and is function of Rp and Tg . For the first drying phase this wet radius equals the radius of the wet particle, and for the second drying phase it corresponds to the radius of the wet core (Mortier, Van Daele, et al., 2012). More details on this model reduction can be found in Mortier, Van Daele, et al. (2012). All the work referred to was performed for a single drying granule. However, the behaviour of a batch of granules is of particular interest since individual granules are likely to influence the ambient conditions, which again influence the drying behaviour of other particles. 2. Population Balance Modelling The development of a complete mechanistic model for the fluidized bed process that allows predicting the spatial changes in the

distribution of the moisture content in a population of granules and eventually the moisture content distribution in the outlet product of the dryer was approached in a step-wise manner. The model for single pharmaceutical granules should first be extended towards a population of granules using the well-known (in engineering, but not in pharmaceutical sciences) framework of Population Balance Model(ling) (PBM). A stand-alone PBM-model is a first step towards a potentially coupled Computational Fluid Dynamics (CFD)-PBMmodel, should this be required. The behaviour of many systems, including a pharmaceutical drying system is often described by the ‘averaged’ behaviour of single particles. However, in a real system this could turn out to be an overly rough approximation. Spatial and population heterogeneity occur in a real system, implying that the ‘state’ of the particles can be different and also the environment they ‘observe’ is not the same throughout the system. This can potentially have an impact on the system behaviour which can be quite different from the ‘average’ behaviour. It really depends on the goal of the modelling study whether this level of detail is required or not. Since the goal of this work is to develop detailed process knowledge, it is important to explore the impact of these heterogeneities. For a thorough analysis of particles that are interacting with each other and the continuous phase, PBMs can be applied (Mortier et al., 2011). If the ambient condition is identical for all granules, a stand-alone PBM-model gives no added-value for understanding the drying process. However, after gaining knowledge about the flow pattern of granules in a fluidized bed, this information can be coupled with the PBM. Individual granules are not subjected to the same ambient condition, so a distinct drying behaviour can be expected. Moreover, the dryer is fed over a longer period of time implying that not all granules are at the same stage in the drying process. This effect can also be accounted for using a PBM. The change of the number density n(x, r, t) can be described by the general Population Balance Equation (PBE):

∂ ˙ ˙ r, Y, t)n(x, r, t) + ∇ R(x, r, Y, t)n(x, r, t) n(x, r, t) + ∇ X(x, ∂t = h(x, r, Y, t)

(1)

where x, r, t, Y and h are respectively the internal coordinate (i.e. the internal property of the particles that is considered to be distributed), the external (spatial) coordinate, the time, the continuous phase vector and the net birth rate due to discrete processes. The net birth rate can include different phenomena such as nucleation, aggregation, agglomeration, breakage, etc. (Ramkrishna, 2000). ∇ , the nabla operator, represents the vector differential operator. In

Fig. 2. Model reduction strategy in order to reduce the full drying model which is too complex to be incorporated in directly (after Mortier, Van Daele et al. (2012)).

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53

a three-dimensional Cartesian coordinate system, ∇ is defined in terms of partial derivative operators as

∇ = xˆ

∂ ∂ ∂ + yˆ + zˆ ∂x ∂y ∂z

(2)

where xˆ , yˆ and zˆ represent the unit vectors in each direction. X˙ and R˙ are the partial derivatives of the internal and external coordinates respectively (Ramkrishna, 2000). In the system at hand the internal coordinate of interest is a representation of the moisture content of the granule, which corresponds to the wet radius of the particle. It is assumed that the particles do not break up during drying. When breakage would be taken into account, the one dimensionial PBM should be extended towards a two dimensional form in order to also consider particle size as an additional dimension in the vector x. Increasing the number of internal coordinates has the disadvantage that it leads to an increase in mathematical complexity and required computational power. Moreover, it also results in difficulties in case one later on would like to combine the PBM with a CFD-code. Therefore it was not considered at this stage. Since drying of wet granules is a continuous process, so no discrete processes occur, h equals zero. When the ambient air phase is assumed to be constant, meaning there is no spatial variation, the third term of the general PBE (Eq. (1)) can be removed. The resulting reduced, homogeneous PBE is a linear hyperbolic equation in (Rw , t).

∂ ∂ ˙ Rw (Rw , Y )n(Rw , t) = 0 n(Rw , t) + ∂t ∂Rw

(3)

where Rw is the wet radius of the particle. In this equation the growth term Gr is defined as Gr = R˙w (Rw , Y )

(4)

In the literature, population balances are typically found with a positive growth term. However, in this paper a negative growth term occurs, similar to a PBM that models dissolution. The difficulty of the incorporation of a growth term in a PBE is that they lead to hyperbolicity in the resulting equations. Since there was no analytical solution for the developed PBM identified, a numerical method to solve the PBE is required for its solution. Several solution methods are available to solve this type of PBEs: Method of Moments (MOM), Monte Carlo simulation, discretization methods, methods of weighted residuals/orthogonal collocation are some examples (Gunawan, Fusman, & Braatz, 2004). The population density function can span over several orders of magnitudes and the distribution can show very steep changes. Both factors complicate achieving an accurate numerical simulation of the PBM within an acceptable computational time (Qamar, Elsner, Angelov, Warnecke, & Seidel-Morgenstern, 2006).

41

into a PBM since the original model is too complex. The empirical equation for the first drying phase was determined to be: Gr,1 (Rw,nor , Tg ) = A + B · Rw,nor + C · eD · Rw,nor Rw,nor =

(5)

Rw − Rp Rw,0 − Rp

(6)

where A, B, C and D are empirical coefficients, which are dependent on Tg , and Rw,0 is the initial (wet) radius (Mortier, Van Daele, et al., 2012). The behaviour of the second drying phase is described by: 

   Gr,2 (Rw,nor , Tg ) = A · (Rw,nor )B + C  · (1 + D · Rw,nor )E





E 

+ Rf ∗ (A · 0.5B + C  · 1 + D · 0.5

 = Rw,nor



)

(7)

Rw Rp

(8)

with A , B , C , D and Rf empirical coefficients, dependent on Tg (Mortier, Van Daele, et al., 2012). 5. Numerical solution methods 5.1. High Resolution Finite Volume scheme High resolution schemes have already been used for the numerical solution of hyperbolic systems (astrophysical flows, gas dynamics, detonation, etc.) to obtain high accuracy on coarse grids and to resolve sharp discontinuities (Qamar et al., 2006). Qamar et al. started from the high resolution semi-discrete scheme of Koren based on a uniform grid (Koren, 1993). The scheme is discrete in the internal coordinate for the space but continuous in time. A finite volume scheme is applied to Eq. (3), and in order to achieve this, the domain of interest of the internal coordinate is subdivided into N subdomains. Rw,i refers to the cell center of cell i and Rw,i represents the cell width. In this context, the application of the cell centered finite-volume discretization to the PBM generates the following expression (Qamar et al., 2006):



R+

w,i

R− w,i

∂n dRw + ((Gr n)R+ − (Gr n)R− ) = 0 w,i w,i ∂t

(9)

+ − In this semi-discrete equation Rw,i and Rw,i are the cell faces and Gr n are the fluxes. One can now define ni (t) as the average value of the number density in each cell, i.e.

ni (t) =

1 Rw,i



R+

w,i

n(Rw , t)dRw

R−

(10)

w,i

3. Objectives

Using this information Eq. (9) can be rewritten as

The objective of this paper is to develop a PBM for the drying of a batch of pharmaceutical granules. Furthermore, different solution methods are tested. The influence of the grid size, which is important for discretization methods, is investigated. Different solution methods are compared in terms of accuracy and calculation time: the High Resolution Finite Volume (HRFV)-scheme, the Method of Characteristics (MOC) and the Quadrature Method of Moments (QMOM).

∂ni 1 ((Gr n)R+ − (Gr n)R− ) = 0, ∀i = 1, 2 . . . N + R w,i w,i ∂t w,i

The computation of the fluxes at the cell-faces determines the accuracy of the scheme. Using an upwind scheme, the cell fluxes can be calculated, and a high order accuracy can be obtained by piecewise polynomial interpolation: (Gr n)R+ = Gr,R+ w,i

4. Reduced drying model The reduced model for the drying of a single granule at different gas temperatures was described in detail by Mortier, Van Daele, et al. (2012). It was developed for the purpose of integrating it

(11)

,  ∈ [−1, 1]

w,i



ni+1 +



1+ 1− (ni − ni+1 ) + (ni+1 − ni+2 ) 4 4

(12)

When  = −1, the result is a second-order accurate fully one-sided upwind scheme. For  = 1 the standard second-order accurate central scheme is obtained. This piecewise polynomial interpolation is

42

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53

described by van Leer (1985). Qamar et al. have chosen a value of 1/3 for : (Gr n)R+ = Gr,R+ w,i



w,i

1 ni+1 + 2

with

ri+

ri+ =

ni − ni+1 +  ni+1 − ni+2 + 

1





2 + ri+ (ni+1 − ni+2 ) 3 3

Rw (13)

(14)

where  is a parameter with a small value to avoid division by zero. To suppress and avoid wiggles and negative solution values the Sweby’s monotonicity theory is used (Sweby, 1984), which leads to the following: (Gr n)R+ = Gr,R+ w,i

w,i



1 ni+1 + (ri+ )(ni+1 − ni+2 ) 2



KO (r) = max 0, min 2r, min

1

VL (r) =

|r| + r |r| + 1

(17) (18)

R,s (r) = max(0, min(1, 2r), min(2, r))

(19)

Discontinuous solutions have second-order accuracy in smooth regions, while first-order accuracy is used in the district of shock discontinuities. The boundary conditions are not straightforward. At the outflow boundary i = 0 (i.e. completely dried granules) the fully one sided upwind scheme is used ( = −1). At the inflow boundary central interpolation is used ( = 1) for i = N − 1. The resulting boundary conditions are: (Gr n) 1 = Gr, 1 2

2

n1 +



1 (n1 − n2 ) , i = 0, 2

xi− 1 (t), xi+ 1 (t) and the growth term Gr (t, Rw ) is substituted by: 2

dRw Gr (t, Rw ) = dt



= i (t)

2

2

n˜ i :=

w,i

w,i



ni+1 +



1 (ri+ (ni − ni−1 ) , i = 1, 2 . . . N − 2 2 (20)

The MOC starts with the discretization of the domain. In this work a regular grid is used, but also irregular or geometric grids can be applied (Qamar & Warnecke, 2007). N is the number of subdivisions in the domain [Rw,min , Rw,max ]. 2

and

xi+ 1 = xmin + ixi

(21)

2

with xi =

xi− 1 + xi+ 1 2

2

2

R

w,i− 1 2

(t)

(27)



Rw,i (t)

˜ w , t)dRw n(R

(28)

i (t)

  Gr,i+ 1 n˜ i d 2 (Rw,i+ 1 (t) − Rw,i− 1 (t))n˜ i = Rw,i (t) dt Rw,i (t) 2 2

(29)

On the left-hand side of Eq. (29) the product rule is applied, and as a result: Gr,i+ 1 n˜ i dn˜ i n˜ i 2 = − (Gr,i+ 1 − Gr,i− 1 ) Rw,i (t) dt 2 2 Rw,i (t)

(30)

The resulting set of Ordinary Differential Equations (ODEs) consists of Eq. (30) combined with: = Gr,i+ 1 , ∀i = 1, 2 . . . N

(31)

2

5.3. Quadrature Method of Moments

5.2. Method of Characteristics

2

dt

 Rw,i+ 1 (t) 2

˜ w , t) n(R

The resulting ODEs can be solved by a standard ODE solver.

The high resolution scheme with  equal to −1 is also used to investigate the differences of both schemes. For both schemes different flux limiting functions are tested.

x 1 = xmin , xN+ 1 = xmax

w

Gr (t, Rw ) ˜ w , t)dRw n(R Rw

1

dt

(Gr n)N = Gin nin , i = N, (Gr n)R+ = Gr,R+

 dR

The Leibniz formula is applied on Eq. (27) and combined with Eq. (28).

2

(Gr n)N− 1 = Gr,N− 1

(26)

∂ ˜ w , t)dRw + n(R ∂t

˜ dR w,i+ 1

1 (nN + nN−1 ), i = N − 1, 2

(25)

over the control volume i (t) = Eq. (24) is integrated 

i (t)

R,m (r) = max(0, min(r, 1))



∂ ∂ Gr (t, Rw ) ˜ w , t) − ˜ w , t) = 0 ˜ w , t) + Gr (t, Rw )n(R n(R n(R Rw ∂t ∂Rw

(16)

A series of flux limiting functions  are available (LeVeque, 1996; Sweby, 1984): van Leer (1974) (Eq. (17)), minmod (Eq. (18)) (Roe, 1983), superbee (Eq. (19)) (Roe, 1983), Chakravarthy and Osher (1983), monotonized central (van Leer, 1977), UMIST (Lien & Leschzines, 1994). . .

(24)

The combination of this definition with the product rule implies:





2 + r, 2 3 3

(23)

˜ Rw ) := Rw n(t, Rw ) n(t,

2

(15)

with  the flux limiting function (Koren, 1993), defined as:



∂ ∂ ˙ Rw (Rw , Y )n(Rw , t) = 0 n(Rw , t) + Rw ∂t ∂Rw

A new variable is defined:

the upwind ratio of two consecutive solution gradients:



The simplified PBE (Eq. (3)) is multiplied by Rw :

and xi = xi+ 1 − xi− 1 2

2

(22)

The method of moments is a solution method where the distribution is approximated by its moments. The kth moment is defined as:



∞ k Rw n(Rw , t)dRw

mk (t) =

(32)

0

The first integral moments are interesting because of their relation to important physical properties. m0 is a representation for the total particle number density (m0 = Nt ), m1 is a kind of average wet diameter (m1 ∼ dt ), m2 represents the total particle wet area (m2 ∼ At ), and m3 the total wet volume (m3 ∼ Vt ). The PBE can be rewritten in terms of moments as: dmk (t) = dt



∞

k−1

kRw Gr (Rw )n(Rw , t)dRw

(33)

0

The moment equations are only closed under certain conditions, which means that the differential equations for the lower order moments do not depend on values of higher order moments (Hulbert & Katz, 1964; Mortier et al., 2011). The QMOM is a variation where the closure problem is solved using a Gaussian quadrature

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53

43

Table 1 Variable types in the reduced drying model and the PBE. Model

Variable types

Status

Symbol

Number

Reduced model

Known

Fixed by problem Fixed by model

Rp , Tg , Rw,0 Part 1: a1–5 , b1–4 , c1–6 , d1–3     Part 2: a1–2 , b1–2 , c1–2 , d1–2 , e1–2 , Rf,1–3

3 18

PBE

Unknown Known

Algebraic Fixed by problem Fixed by model

Gr,1 , Gr,2 Rp , Tg , , Part 1: a1–5 , b1–4 , c1–6 , d1–3     , d1–2 , e1–2 , Rf,1–3 Part 2: a1–2 , b1–2 , c1–2

Unknown

Differential

Rw , n

approximation (Marchisio, Pikturna, Fox, Vigil, & Barresi, 2003; McGraw, 1997): n(Rw , t) ≈

Nq

wi (t)[Rw − Rw,i (t)]

(34)

i=1

Using this quadrature approximation the moments become: mk (t) ≈

Nq

k wi (t)Rw,i (t)

(35)

i=1

where Rw,i (t) are the abscissas and wi (t) are the weights, which can be specified from the lower order moments. The weights and abscissas can be determined by different algorithms; the Product-Difference (PD)-algorithm and the Chebyshev algorithm. A sequence of moments is used by the PD-algorithm for the computation of the coefficients of the ‘three-term-recurrence relationship’ that is satisfied by orthogonal polynomials. With the help of these coefficients the quadrature points and weights are calculated, which are then used to approximate the integrals over the unknown number density function. The Chebyshev algorithm had been used previously for the computation of the coefficients of the recurrence relations of orthogonal polynomials with known density function. Upadhyay adapted the algorithm to work with a sequence of moments with the underlying density function being unknown. The algorithm uses directly the ‘three-term-recurrence relationship’ among the orthogonal polynomials to establish recursive formulae for the coefficients (Upadhyay, 2012). According to Upadhyay the Chebyshev algorithm is more general and more robust, and it can also be applied to symmetric distributions, where the PD-algorithm would struggle with odd moments which have a value of zero. High order moments can be calculated to obtain accurate high order quadratures (Upadhyay, 2012). 6. Model analysis and simulation parameters In Table 1 the variable types are listed for the reduced drying model, in which a division is made between the variables which should be specified (the known variables) and the variables predicted with the model equations (unknown variables). For the reduced model the growth term of the first and second drying phase are the unknowns, and are calculated using algebraic equations (Eqs. (5)–(7)). The reduced drying model consists of respectively 18 and 13 coefficients for the first and second drying phase. The gas temperature, particle radius and initial moisture content are fixed in the problem to be solved (Mortier, Van Daele, et al., 2012). Simulations for testing the PBM were performed at a specific gas temperature and for particles with a certain size (Table 2). A normal probability density function was used as initial distribution of the internal coordinate: n0 =

2 1 2 √ e−(Rw −) /(2 ) 2

(36)

Total

13 2 4 18

34 2

13 2

35 2

with n0 the initial distribution, the standard deviation and  the mean. The known and unknown variables of the used PBE are listed in Table 1. As the reduced model is used to calculate the growth term of the PBE the number of coefficients is the same for both models. Using the PBE two variables are calculated Rw , n. For the HRFV-scheme the discretization domain ranged from 0 till 1.032 times Rp . The MOC uses a moving domain, so the minimum initial value of the domain should not be 0, but was set to Rp . The maximum value for the initial discretization domain was also set at 1.032 times Rp . All calculations were performed using Matlab®, and a build-in ODE-solver was used (ode23, which uses a variable step size for the time). 7. Results 7.1. Solution methods based on discretization of the domain The two methods based on discretization are first evaluated for a very small grid size, which is assumed to represent the pseudoanalytical solution. For the HRFV-scheme this is done for a -value equal to 1/3 and −1 (Sections 7.1.1, 7.1.4 and 7.1.6). The influence of decreasing the grid size is investigated for the HRFV-scheme with a -value of 1/3 and the MOC (Sections 7.1.2 and 7.1.7). Furthermore the effect of different flux limiting functions is analyzed for the HRFV-scheme (Sections 7.1.3 and 7.1.5). In Section 7.1.8 the different methods based on discretization are compared and conclusions are drawn. 7.1.1. PBM solution using the HRFV-scheme with highest accuracy for a -value of 1/3) The HRFV-scheme was in a first approach used with the largest accuracy (i.e. N = 2000). This is assumed to be the most accurate solution that can be obtained within a reasonable calculation time. A too high grid size (e.g. 5000) lead to a significantly larger computational burden (calculation time of several weeks). Moreover, further increasing N (>2000) did not result in an improvement of the accuracy. In Figs. 3–5 the model predictions are presented for N equal to 2000 and the flux limiting function of Koren (KO ). The normalised moments are calculated using the discrete version of Eq. (32) (normalisation based on moment value at t = 0). At the beginning the curve of m0 displays a dynamic behaviour, but after less than 5 s, the curve becomes horizontal (Fig. 3). As the total number of particles does not change in time, m0 should ideally not change in time and Table 2 Simulation parameters. Rp Tg 

6e−4 m 55 ◦ C 1.02 Rp 3e−6

44

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 1.0015

4

14

x 10

t=5s t=7s t=9s t = 11 s

12

10

1.001 1.002

8

1.0005

n

m0

m0

1.0015 1.001 1.0005 1

6

4 0

5

10 t (s)

15

20

2

0 1

0

50

100

150

200 t (s)

250

300

350

400

−2

5.6

5.7

5.8

Fig. 3. m0 (t) for the HRFV-scheme with  = 1/3 and N = 2000.

0.9

1

0.8 m1

0.95

0.7

m1

6

6.1

6.2 −4

x 10

Fig. 6. n(Rw ) at different time instants for the HRFV-scheme with  = 1/3 and N = 2000 with the focus on the bimodal distribution.

1

0.9

0.6

0.85

0

5

10 t (s)

15

20

0.5

0.4

0.3

0.2

5.9 Rw (m)

0

50

100

150

200 t (s)

250

300

350

400

Fig. 4. m1 (t) for the HRFV-scheme with  = 1/3 and N = 2000.

stay at 1. This is established apart from a numerical error of 0.13%. The first moment, a representation of the average wet diameter, should logically decrease from the start which is indeed observed (Fig. 4). In Fig. 5 the number density is presented at different points in time. It is clear that the wet radius decreases. Because some particles complete the first drying phase sooner compared to other particles, the width of the number density distribution increases a bit as time evolves. However, this effect is limited as all particles are subjected to the same ambient conditions. In Fig. 6 the transient between drying phase 1 and 2 is illustrated in more detail for a population of particles. It can be seen that this effect explains the transient observed in the inset of Fig. 4 and only acts during the first few seconds of the simulation (depending on Tg ). In Fig. 7 the standard deviation of the moisture distribution is plotted as a function of time. A clear increase in the beginning can be seen, which represents the occurrence of a bimodal distribution. This is followed by a decrease (bimodal distribution again becomes a monomodel distribution), reaching a minimum value. Later on, −6

11

x 10

10 9

Standard deviation

8 7 6 5 4 3 2 1

Fig. 5. n(Rw ) at different time instants for the HRFV-scheme with  = 1/3 and N = 2000 (the value for Rw decreases in time).

0

50

100

150

200 t (s)

250

300

350

400

Fig. 7. The standard deviation of the distribution as a function of time for the HRFVscheme with  = 1/3 and N = 2000.

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 1

1000

45

1.025

Mean RMSE tCal

KO VL Rm Rs

1.021

1.015 500

1.02

m0

0.5

tCal

meanRMSE

1.02

1.019 1.01

1.018 1.017 1.016

1.005

0 50

100

150

200

250

300

350

400

450

0 500

N

Fig. 8. The RMSE and the calculation time for increasing values of N for the HRFVscheme with  = 1/3.

an increase can again be observed which means the variance of the distribution widens again due to the difference in growth rate of particles with different moisture content (cfr. size dependent growth term). At the end of the drying process the wet radius is very small, and for a certain decrease in wet radius only a limited amount of water should evaporate. This is why at the end of the simulation the evaporation rate increases again, but as some particles still have a larger wet radius, with a lower evaporation rate, the width of the distribution increases at the end. 7.1.2. Investigation of grid coarseness for the HRFV-scheme ( = 1/3) To investigate the effect of the grid size N, as well as to determine the ‘optimal’ grid size, N was decreased and varied between 50 and 500. The scenario with an N-value of 2000 was used as reference for the calculation of the Root Mean Squared Error (RMSE) for the other grid sizes assuming that to be the pseudo-analytical solution. The RMSE was calculated as the mean of the first three moments. This was done using KO as flux limiting function. The effect on the calculation time of the decrease in N was also investigated. In Fig. 8 the RMSE and the calculation time are visualised. Starting from a value of 300 for N the curve of the RMSE flattens and there is almost no improvement. In Table 3 the calculation time is listed for all simulations with different grid sizes, where the calculation time is normalised with respect to the smallest obtained value (i.e. N = 50). Also the calculation time for a grid size of 2000 is mentioned. The calculation time is scaled to allow better comparison of calculation times obtained with different solution techniques. It is obvious that the calculation time increases with N. The calculation time is however not linearly Table 3 Scaled calculation time for the HRFV-scheme using different N-values for  = 1/3. N 2000 500 450 400 350 300 250 200 150 100 50

Scaled calculation time 1.21e4 190 147 112 79.1 56.4 34.7 20.3 10.71 4.27 1

1

0

50

100

150

200 t (s)

2

3

250

4

300

5

6

350

400

Fig. 9. m1 (t) for the HRFV-scheme with  = 1/3 and N = 300 with different flux limiting functions.

dependent on N. Comparing a grid size of 2000 with a grid size of 200 gives a reduction of the calculation time with more than 99.8%. At N = 300, a reduction of 99.5% is obtained. Comparing the effect of different grid sizes on the moments and the number density for a -value of 1/3, the following conclusions can be made: As N increases, the height of the peak at the start of m0 decreases (Fig. S1 (Supplementary data)). During the first stages of the drying process, a significant difference can be seen between the moment predictions, especially for coarse grids (N < 250). The time needed for m0 to reach steady state gradually increases for lower N. In Fig. S2 (Supplementary data) the first moment is presented for different grid sizes. For N < 250, significantly deviating behaviour is observed, whereas, the other grid sizes display more or less the same trend. In Fig. S3 (Supplementary data) the number density is presented after 300 s. The width of the distribution increases clearly with decreasing N-value. This is also true for N > 250. However, these deviations do not seem to significantly impact the lower integral moments of the distribution. Depending on what one’s interest is, one needs to be careful in choosing N. A value larger than 2000 might even be required. The standard deviation, as a measurement of the width of the distribution decreases clearly with increasing values of N (Fig. S4 (Supplementary data)). For small values of N (< 200) the curve of the standard deviation shows some sharp peaks at later time steps. 7.1.3. Influence of the flux limiting function for the HRFV-scheme ( = 1/3) The influence of the flux limiting function is investigated in more detail for an N-value of 300, as higher grid sizes did not improve the results significantly. In Fig. 9 the 0th moment is presented, and for the early time steps a difference can be detected. The number density is also affected by this and at 300 s the width is the smallest when using R,s , followed by KO , VL and R,m (Fig. 10). In Fig. 11 the standard deviation is presented in function of time for the different flux limiting functions. The difference in width of the distribution presented in Fig. 10 is clearly visible. The width of the distribution using R,s as flux limiting function is obviously smaller compared to the others. In Table 4 the scaled calculation time is listed. The sequence of the calculation time is the same as that of the standard deviation (i.e. shorter calculation times coincide with larger final standard deviation).

46

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 4

8

1.009

x 10

KO VL Rm Rs

7

1.008 1.007

6 1.006

5

1.01 1.005

1.008

m0

4

m0

n

1.004

3

1.006 1.004

1.003

1.002

2 1

1.002

0

5 t (s)

1

10

1.001

0 1

−1

2.1

2.2

2.3

2.4

2.5 Rw (m)

2.6

2.7

2.8

2.9

−6

x 10

KO VL Rm Rs

14

Standard deviation

12

10

8

6

4

2

0

50

100

150

200 t (s)

250

300

350

50

100

150

200 t (s)

250

300

350

400

3 −4

Fig. 12. m0 (t) for the HRFV-scheme with  = −1 and N = 2000.

x 10

Fig. 10. n(Rw ) at 300 s for the HRFV-scheme with  = 1/3 and N = 300 with different flux limiting functions.

16

0

400

steps showed negative values, which is physically impossible and is caused by numerical instability (Fig. 14). This renders these simulations highly uncertain. The standard deviation is presented in Fig. 15. Near the end of the simulation, values of the standard deviation keep decreasing, meaning that the width of the distribution becomes smaller and smaller. This phenomenon is likely caused by the occurrence of negative values for the number density and is likely a numerical artefact. Investigating the impact of lowering the grid size indicated that the oscillation in the number density at a certain time point was lower for small values of N. However, the negative values for the number density did not dissappear. Due to the uncertainty of these simulations, this method was not further studied. 7.1.5. Influence of the flux limiting function for the HRFV-scheme ( = −1) Again different flux limiting functions are tested for an N-value of 300. The results are presented in Figs. 16–18. It is obvious that the negative values for the number density dissappear when a flux limiting function is implemented. The difference between the methods is limited for m0 . Again the width of the distribution after

Fig. 11. Standard deviation for the HRFV-scheme with  = 1/3 and N = 300 with different flux limiting functions.

1.2 1.1 1.02 1 1 m1

0.9 0.8 m1

7.1.4. PBM solution using the HRFV-scheme with highest accuracy for a -value of −1) An identical evaluation (N = 2000) was performed for a  value of −1 without any flux limiting function. In this case, the result is similar (Fig. 12). Again there is some dynamic behaviour at the start for m0 (Fig. 12). For m1 an initial increase can be observed, which is physically not possible as the moisture content should decrease continuously in time (Fig. 13). The number density at different time

0.98 0.96

0.7

0.94

0

5 t (s)

0.6

Table 4 Scaled calculation time for the HRFV-scheme using different flux limiting functions (N = 300). N KO VL R,m R,s

10

0.5 0.4

Scaled calculation time  = 1/3

 = −1

0.3

4.00 2.46 1 4.76

3.94 2.49 1 4.82

0.2

0

50

100

150

200 t (s)

250

300

350

Fig. 13. m1 (t) for the HRFV-scheme with  = −1 and N = 2000.

400

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53

47

1.025 KO VL Rm Rs

1.02

1.021

1.015 m0

1.02 1.019 1.01

1.018 1.017 1.016

1.005

1

Fig. 14. n(Rw ) at different time instants for the HRFV-scheme with  = −1 and N = 2000 with the focus on the bimodal distribution (the value for Rw decreases in time).

300 s is dependent on the flux limiting function. Comparing the calculation time of the different flux limiting functions, mentioned in Table 4, it is clear that the sequence of the calculation time is the same as that of the standard deviation. In Fig. 18 the standard deviation is presented for the different scenarios. The influence of the flux limiting function on the width of the distribution is significant. 7.1.6. PBM solution using MOC with highest accuracy Figs. 19–21 show the result for the MOC-method using N = 1000. For the normalised m0 a clear peak can be seen at the start, and as time evolves the curve shows a further dynamic behaviour, which becomes more significant at the end (Fig. 19). In fact the value for m0 should stay fixed at 1. Deviations mean that mass leaks have occured. m1 decreases over the whole time range, which would be expected as the moisture content decreases (Fig. 20). The width of the number density at different time steps is only slightly different, but does not change much over time (and as such also not the heigth

0

50

100

150

200 t (s)

2

3

4

250

5

300

6

350

400

Fig. 16. m1 (t) for the HRFV-scheme with  = −1 and N = 300 with different flux limiting functions.

of the peak) (Fig. 21). In Fig. 22 the focus is made on the bimodal behaviour at the transition of the first drying phase to the second drying phase. This effect is important to explain the increase of the standard deviation at the start of the drying process. In Fig. 23 the evolution of the standard deviation is presented as a function of time. At the end the increase in width is limited. 7.1.7. Investigation of grid coarseness for the MOC Also for MOC the effect of N was investigated. Due to the moving domain the grid size can be much lower compared to the HRFVscheme. In this case the grid size comprised values ranging from 10 to 1000 (Table 5). The decrease in RMSE with increasing N is presented in Fig. 24, and a significant increase in calculation time can be observed. The decrease in RMSE is calculated using N equal to 1000 as the reference case. The scaled calculation time is mentioned in Table 5. It is obvious that the calculation time increases with N, but no linear dependency can be observed.

4

−6

11

8

x 10

KO VL Rm Rs

7

10 9

6

8

5

7

4 n

Standard deviation

x 10

6

3

5

2 4

1 3

0 2 1

−1 0

50

100

150

200 t (s)

250

300

350

400

Fig. 15. The standard deviation of the distribution as a function of time for the HRFV-scheme with  = −1 and N = 2000.

2.1

2.2

2.3

2.4

2.5 Rw (m)

2.6

2.7

2.8

2.9

3 −4

x 10

Fig. 17. n(Rw ) at 300 s for the HRFV-scheme with  = −1 and N = 300 with different flux limiting functions.

48

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 −6

x 10

16

KO VL Rm Rs

14

Standard deviation

12

10

8

6

4

2

0

50

100

150

200 t (s)

250

300

350

400

Fig. 18. Standard deviation for the HRFV-scheme with  = −1 and N = 300 with different flux limiting functions.

Fig. 21. n(Rw ) at different time instants for MOC with N = 1000 (the value for Rw decreases in time).

1.005 1.004 4

1.003

14

1.002

12

x 10

t=5s t=7s t=9s t = 11 s

m0

1.001

10 1 0.999

n

8

0.998

6 0.997

4 0.996

0

50

100

150

200 t (s)

250

300

350

400

2

Fig. 19. m0 (t) for the MOC with N = 1000. 0

1

5.6

5.7

5.8

5.9 Rw (m)

0.9

0.7

Table 5 Scaled calculation time for MOC using different N-values. m1

6.1

6.2 −4

x 10

Fig. 22. n(Rw ) at different time instants for MOC with N = 1000 with the focus on the bimodal distribution.

0.8

0.6

N 0.5

0.4

0.3

0.2

6

0

50

100

150

200 t (s)

250

300

Fig. 20. m1 (t) for the MOC with N = 1000.

350

400

1000 500 400 300 200 100 50 40 30 20 10

Scaled calculation time 6.34e3 1.54e3 1.06e3 594 263 68.1 18.1 11.76 7.02 3.36 1

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 −5

1.2

49

1.045

x 10

HRFV 1/3 KO HRFV 1/3 Rs MOC

1.04 1.035

1

1.03 1.03

1.025 m0

Standard deviation

0.8

1.02

1.02

0.6

1.015 1.01 1.01

0.4

1

1.005

0

1

2

3

1

0.2

0.995 0 0

50

100

150

200 t (s)

250

300

350

400

0

50

100

150

200 t (s)

250

300

350

400

Fig. 25. Comparison of m0 (t) for the different solution methods using discretization.

Fig. 23. The standard deviation of the distribution as a function of time for the MOC.

For m0 a clear peak can be seen at the start for each grid size (Fig. S5 (Supplementary data)). This peak does not disappear for larger N, but decreases significantly. At the end a constant difference can be seen between the grid sizes. For m1 only a clear difference can be detected at the start when comparing different grid sizes (Fig. S6 (Supplementary data)). The number density after 300 s is almost the same for all grid sizes, but the curve for N equal to 10 is less smooth (Fig. S7 (Supplementary data)). The overall model behaviour seems to remain unchanged as of N = 100, which seems the optimal value for this solution method (trade off between the required calculation time and the obtained accuracy). 7.1.8. Comparison and conclusion for the methods based on discretization of the domain It is interesting to compare the HRFV-schemes mutually and with the MOC in terms of the number density, the width of the distribution, the calculation time and the moment dynamics. Without implementing a flux limiting function in the HRFVscheme the result is significantly different (occurrence of negative values for the number density when  equal to −1). However, after implementation of a flux limiter in both HRFV-schemes still some differences can be observed. The sequence of the standard deviation is comparable for both -values (Figs. 11 and 18). Comparing 0.01

1000

0.005

0

500

0

50

100

150

200

250 N

300

350

400

450

tCal

meanRMSE

Mean RMSE tCal

0 500

Fig. 24. The RMSE and the calculation time for increasing values of N for the MOC.

the absolute values for the width of the distribution it is obvious that the result using VL and KO are closer to each other for a value of −1 compared to a -value of 1/3. For VL , R,m and R,s the width is almost identical for both -values. The ranking of the calculation times is the same for both -values (Table 4). The number density after 300 s is, similar to the standard deviation, only significantly different between both -values for KO (Figs. 10 and 17). Next, the HRFV-scheme using KO (as proposed by Qamar et al., 2006) and R,s (which results in the smallest distribution) can be compared with MOC. This was done using a grid size of 300 for the HRFV-scheme, as there is almost no improvement in the RMSE when a higher grid size is used (Fig. 8). For MOC N was set equal to 100, since a higher grid size has almost no influence on the number density at 300 s (Fig. S7 (Supplementary data)). In Fig. 25 the zeroth moment is presented. The result for the HRFV-scheme is somewhat different compared to MOC. At the start both methods gave an increase followed by a decrease. However, the result for the HRFV-scheme remains stable after the initial dynamic behaviour whereas the result for the MOC fluctuates also at later time points. The m1 should decrease from the start, which happened for MOC, but this is not the case for the HRFV-scheme (Fig. 26). The results for the higher order moments are presented in Figs. S8–S11 (Supplementary data). For m2 till m5 a small bump can be observed in the beginning for the HRFV-scheme, while this is not the case for the MOC. In Fig. 27 the initial distribution is presented. The results for both flux limiting functions used in the HRFV-scheme gave of course the same result, but as the discretization of the domain is different for the MOC a small difference could be expected. At 120 s a clear difference in the width and height of the peak can already be detected (Fig. 28). Where the solution is a sharp peak for the MOC, the HRFVschemes give a curve which is much wider and exhibits a lower peak. At higher time steps (320 and 400 s) the distributions resulting from the HRFV-scheme become wider and wider whereas that of the MOC almost remains unaltered (Figs. 29 and 30). In Fig. 31 the standard deviation is presented. It is obvious that the width of the distribution is strongly dependent on the solution method. While the increase in the standard deviation is limited for the MOC, the increase for the HRFV-scheme is more pronounced. The distribution at the end is clearly smaller when the MOC is used. For the HRFV-scheme the width of the distribution is strongly dependent on the grid size (Fig. S4 (Supplementary data)). An almost identical curve for the standard deviation could be obtained

50

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 5

1.2

2.5

HRFV 1/3 KO HRFV 1/3 Rs MOC

1.1 1

x 10

HRFV 1/3 KO HRFV 1/3 Rs MOC

2

0.9

5

0.8

1.03

0.7

1.02

0.6

1.01

1.5

x 10

3

n

m1

2 1 0 1

0.5

0.5 −1 0.99

0.4

0

1

2

3.6

3.8

4

4.2

3

4.4 −4

x 10 0

0.3 0.2

1

0

50

100

150

200 t (s)

250

300

350

400

−0.5

1.5

2

2.5

3

3.5 4 Rw (m)

4.5

5

5.5

6

6.5 −4

x 10

Fig. 26. Comparison of m1 (t) for the different solution methods using discretization.

Fig. 28. n(Rw ) at t = 120 s for the HRFV-scheme with  = 1/3 (KO and R,s ) and the MOC-method.

when N would be chosen around 2000. However, this requires significantly longer calculation times. But based on this fact, it could be concluded that the MOC approximates the number density more precisely, as an identical result would have been obtained when the accuracy of the HRFV-scheme would have been increased (by increasing N). The calculation time is definitely different for all methods. Whereas only 21 s were required for the MOC, the HRFV-schemes took 284 and 277 s for KO and R,s respectively. Based on the performance of the distinct discretization methods and the calculation time, the MOC seems to be the preferred method for the PBM at hand. The HRFV-scheme with a -value of −1 without flux limiting function leads to negative values for the number density, which is not desirable. Increasing N did not give better results than low values for N, and in fact the oscillation was even lower for low values of N.

7.2. Quadrature Method of Moments The QMOM is investigated with the PD- and the Chebyshev algorithm. The same simulation parameters were used as for the discretization methods, and the grid size to calculate the initial distribution in terms of n was set to 500. For this method a large N for the calculation of the initial distribution gives more accurate values for the initial moments, but has no influence on the total calculation time. For the other solution methods the grid size for the initial distribution is chosen according to the grid size used during calculation. In this case the number of moments that are calculated is varied between 3 and 7. The results for the PD-algorithm are presented in Figs. 32–34 for the first three moments and in Figs. S12–S15 (Supplementary data) for the higher order moments. With the PD-algorithm it was possible to calculate 7 moments. The difference between the higher moments when 3 or more moments are calculated is minimal.

4

14

x 10

5

3

HRFV 1/3 KO HRFV 1/3 Rs MOC

12

x 10

HRFV 1/3 KO HRFV 1/3 Rs MOC

2.5

5

10

2

3

x 10

4

15

x 10

2

1.5

1

n

n

8

10 6

1

0

5 −1

4

0.5

6

6.2

2.6 x 10

−4

x 10

1.5

2.4

6.4

2

0

2.2

−4

0

2

2.5

3

3.5 4 Rw (m)

0

4.5

5

5.5

6

6.5 −4

x 10

Fig. 27. n(Rw ) at t = 0 s (initial condition) for the HRFV-scheme with  = 1/3 (KO and R,s ) and the MOC-method.

−0.5

1.5

2

2.5

3

3.5 4 Rw (m)

4.5

5

5.5

6

6.5 −4

x 10

Fig. 29. n(Rw ) at t = 320 s for the HRFV-scheme with  = 1/3 (KO and R,s ) and the MOC-method.

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 5

2.5

51

2

x 10

HRFV 1/3 KO HRFV 1/3 Rs MOC

2

3 4 5 6 7

1.5 5

3

1.5

x 10

1 m0

2 n

1 1

2 0.5 1

0.5

0

1.4

1.6

1.8 −4

0

0

x 10 0

−1 −0.5

−0.5

1.5

2

2.5

3

3.5 4 Rw (m)

4.5

5

5.5

6

0

5

0

50

100

10 150

200 t (s)

6.5

250

300

350

400

−4

x 10

Fig. 32. m0 (t) for the QMOM using the PD-algorithm. Fig. 30. n(Rw ) at t = 400 s for the HRFV-scheme with  = 1/3 (KO and R,s ) and the MOC-method.

−4

7

x 10

3 4 5 6 7

−5

x 10

1.2

6

HRFV 1/3 KO HRFV 1/3 Rs MOC

1

0.8

m1

Standard deviation

5

−4

4

6.1

x 10

6.05

0.6

3 6 0.4

5.95

2

5.9 0.2

1 0

0

50

100

150

200 t (s)

250

300

350

2

0

50

400

100

6 150

200 t (s)

nMom

Scaled calculation time

3 4 5 6 7

0.33 0.49 0.65 0.88 1

300

350

400

−7

4

x 10

3 4 5 6 7

3.5 3 −7

3.8

m2

2.5

x 10

3.7 2 3.6 1.5 3.5

Table 6 Scaled calculation time for QMOM using the PD-algorithm.

250

Fig. 33. m1 (t) for the QMOM using the PD-algorithm.

Fig. 31. Comparison of the standard deviation for the different solution methods using discretization.

Also the increase in calculation time is limited when increasing the number of moments (Table 6). Using the Chebyshev algorithm the number of moments was limited to 3 due to infinity problems. The moments are presented in Figs. S16–S18 (Supplementary data). Upadhyay concluded that the Chebyshev algorithm can be used for very high order moments to obtain accurate high order quadratures (Upadhyay, 2012). However, based on the presented results the PD-algorithm was able to calculate at least 7 moments, whereas the Chebyshev algorithm

4

1

2

4

6

0.5 0

0

50

100

150

200 t (s)

250

300

Fig. 34. m2 for the QMOM using the PD-algorithm.

350

400

52

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53 Table 7 Scaled calculation time to compare a discretization method with a method using moments.

1.025 MOC QMOM PD 1.02

1.015

Method

Scaled calculation time

MOC QMOM-PD

45.55 1

m0

1.025 1.02

1.01

1.015 1.01 1.005

1.005 1

0

2

4

6

1

0.995

0

50

100

150

200 t (s)

250

300

350

400

8. General conclusion

Fig. 35. Comparison of m0 (t) for the MOC-method with QMOM-PD.

failed after more than 3 moments. The calculation time when using the Chebyshev algorithm was 6 times the calculation time when using the PD-algorithm. 7.3. Comparing MOC (discretization) and QMOM-PD (moments-based method) Two different types of solution methodologies are compared, one is based on discretization (MOC), while the other is based on moments (QMOM-PD). For the MOC N was set to 100. To compare the discretization methods with the method of moments the moments were calculated with the number density at each time step. The zeroth moment shows a difference between the discretization method and the method of moments (Fig. 35). The MOC results in a sharp peak at the start and a bump at the end. For m1 a difference is observed at the start, but the trend for both methods is the same (Fig. 36). The detailed part of the figure shows that the method of moments develops approximately the same. The same conclusion can be drawn for m2 (Fig. S19 (Supplementary data)). The other moments give approximately the same result (Figs. S20–S22 (Supplementary data)). In fact it can be concluded

1 MOC QMOM PD 0.9

0.8

m1

0.7

1

0.6

0.99

0.5

0.98

0.4

0.97

0

2

4

6

0.3

0.2

0

50

100

150

200 t (s)

250

300

350

Fig. 36. Comparison of m1 (t) for the MOC-method with QMOM-PD.

that both methods would end up with the same results in terms of the moments. The calculation time for the methods of moments is significantly lower (Table 7), but the result is a set of moments compared to the full number density for the discretization methods. It is important to mention that QMOM-PD is very fast to run, i.e. it took only half a second to run the simulation. As such, it is an interesting tool to extend with other models, such as CFD-models. In the literature several methods for the reconstruction of the moments are available (John, Angelov, Öncül, & Thévenin, 2007).

400

The implementation of the reduced drying model in a PBM enables to calculate the evolution of the moisture content for a batch of granules. The PBE was solved using different solution methodologies: • The HRFV-scheme with a -value equal to 1/3 calculates the number density accurately using a grid size higher than 300. A grid size of 50 results in a significant difference in the moments. However, the width of the number density decreases when a larger grid size is used. A grid size of 300 was chosen for further research. The choice of the flux limiting function does not have a large influence on the result, only the calculation time can be reduced by approximately a factor of 5 using R,m in stead of R,s . Using a -value equal to −1 without any flux limiting function it was not possible to calculate the distribution correctly, as negative values for the number density occur. This problem can be solved by implementing a flux limiter. • The MOC, which uses a moving grid, requires a lower grid size to obtain a smooth curve for the number density. Here, a grid size of 100 was sufficient. • Comparing the methods using a discretization of the domain it can be concluded that the MOC is much faster than the HRFVscheme with a -value of 1/3. The MOC uses a moving grid, and as such the grid size can be taken lower and still the result will be accurate. • The QMOM with the PD-algorithm is able to calculate at least 7 moments. The increase in calculation time is limited when increasing the number of moments. When using the Chebyshev algorithm it is only possible to calculate 3 moments, which is in contradiction with the result presented by Upadhyay (2012). • QMOM is significantly faster than the discretization methods but no number density is obtained as a result. There exist several reconstruction methods but when using such methods the computation time increases again. The PBM-model with the implemented empirical drying model allows to calculate the evolution of the moisture content during drying. The model can now be used to compute the evolution of the distribution for different initial distributions, different particle sizes and different gas temperatures. The model can be used to investigate the effect of the change in gas temperature during drying, and as such to develop some guidelines about the set-up of the dryer. However, before using the model for process control and defining the design space, a validation should be carried out using experimental data.

S.T.F.C. Mortier et al. / Computers and Chemical Engineering 50 (2013) 39–53

The used modelling approach, which is first setting up and calibrating/validating a drying model for a single particle, followed by model reduction and inclusion in a PBM, can definitely be generalized. However, the model formulation of the PBM might turn out to be different as well as parameter values obtained. Acknowledgements Financial support for this research from the Fund for Scientific Research Flanders (FWO Flanders – Ph.D. fellowship Séverine Thérèse F.C. Mortier) is gratefully acknowledged. Appendix A. Supplementary Data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.compchemeng .2012.11.005. References Chakravarthy, S., & Osher, S. (1983). High resolution applications of the Osher upwind scheme for the Euler equations. In 6th CFD conference (pp. 363–372). Gernaey, K., Cervera-Padrell, A., & Woodley, J. (2012). A perspective on process systems engineering in pharmaceutical process development and innovation. Computers and Chemical Engineering, 42, 15–29. Gunawan, R., Fusman, I., & Braatz, R. (2004). High resolution algorithms for multidimensional population balance equations. AIChE Journal, 50(11), 2738–2749. http://dx.doi.org/10.1002/aic.10228. URL http://doi.wiley.com/10.1002/aic.10228 Hegedus, A., & Pintye-Hódi, K. (2007). Comparison of the effects of different drying techniques on properties of granules and tablets made on a production scale. International Journal of Pharmaceutics, 330(1–2), 99–104. http://dx.doi.org/10.1016/j.ijpharm.2006.09.001 http://www.ncbi.nlm.nih.gov/pubmed/17049769 Hulbert, H., & Katz, S. (1964). Some problems in particle technology: A statistical mechanical formulation. Chemical Engineering Science, 19(8), 555–574. John, V., Angelov, I., Öncül, A., & Thévenin, D. (2007). Techniques for the reconstruction of a distribution from a nite number of its moments. Chemical Engineering Science, 62, 2890–2904. Koren, B. (1993). A robust upwind discretisation method for advection, diffusion and source terms. In Numerical methods for advection-diffusion problems. Friedrich Vieweg & Sohn Verlag. (pp. 117–138). LeVeque, R. (1996). High-resolution conservative algorithms for advection in incompressible flow. SIAM Journal on Numerical Analysis, 33, 627–665. Lien, F., & Leschzines, M. (1994). Upstream monotonic interpolation for scalar transport with application to complex turbulent flows. International Journal for Numerical Methods in Fluids, 19(6), 527–548.

53

Marchisio, D., Pikturna, J., Fox, R., Vigil, R., & Barresi, A. (2003). Quadrature method of moments for population-balance equations. AIChE 49(5), 1266–1276. http://dx.doi.org/10.1002/aic.690490517 Journal, http://doi.wiley.com/10.1002/aic.690490517 McGraw, R. (1997). Description of aerosol dynamics by the quadrature method of moments. Aerosol Science and Technology, 27(2), 255–265. http://dx.doi. org/10.1080/02786829708965471, http://www.informaworld.com/openurl? genre=article&doi=10.1080/02786829708965471&magic=crossref||D404A21C 5BB053405B1A640AFFD44AE3 Mortier, S., De Beer, T., Gernaey, K., Remon, J., Vervaet, C., & Nopens, I. (2011). Mechanistic modelling of Fluidized Bed Drying processes of wet porous granules: A review. European Journal of Pharmaceutics and Biopharmaceutics, 79, 205–225. Mortier, S., De Beer, T., Gernaey, K., Vercruysse, J., Fonteyne, M., Remon, J., et al. (2012). Mechanistic modelling of the drying behaviour of single pharmaceutical granules. European Journal of Pharmaceutics and Biopharmaceutics, 80(3), 682–689. http://dx.doi.org/10.1016/j.ejpb.2011.12.010 Mortier, S., Van Daele, T., De Beer, T., Gernaey, K., Remon, J., Vervaet, C., et al. (2012). Reduction of a single granule drying model: An essential step in preparation of a PBM with a continuous growth term. AIChE Journal, in press, http://dx.doi.org/10.1002/aic.13907 Qamar, S., Elsner, M., Angelov, I., Warnecke, G., & Seidel-Morgenstern, A. (2006). A comparative study of high resolution schemes for solving population balances in crystallization. Computers and Chemical Engineering, 30(6–7), 1119–1131. http://dx.doi.org/10.1016/j.compchemeng.2006.02.012. URL http://linkinghub.elsevier.com/retrieve/pii/S0098135406000445 Qamar, S., & Warnecke, G. (2007). Numerical solution of population balance equations for nucleation, growth and aggregation processes. Computers and Chemical Engineering, 31(12), 1576–1589. http://dx.doi.org/10. 1016/j.compchemeng.2007.01.006 http://linkinghub.elsevier.com/retrieve/pii/ S0098135407000208 Ramkrishna, D. (2000). Population balances: Theory and applications to particulate systems in engineering. New York: Academic Press. http://www. amazon.com/Population-Balances-Applications-Particulate-Engineering/dp/ 0125769709 Roe, P. (1983). Some contributions to the modelling of discontinuous flows. In AMS/SIAM Seminar. Sweby, P. (1984). High resolution schemes using flux limiters for hyperbolic conservation laws. SIAM Journal on Numerical Analysis, 21, 995–1011. Upadhyay, R. (2012). Evaluation of the use of the Chebyshev algorithm with the quadrature method of moments for simulating aerosol dynamics. Journal of Aerosol Science, 44, 11–23. http://dx.doi.org/10.1016/j.jaerosci.2011.09.005 http://linkinghub.elsevier.com/retrieve/pii/S0021850211001613 van Leer, B. (1974). Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second order scheme. Journal of Computational Physics, 14, 361–370. van Leer, B. (1977). Towards the ultimate conservative difference scheme. III. Upstream-centered finite-difference schemes for the ideal compressible flow. Journal of Computational Physics, 23, 263–275. van Leer, B. (1985). Upwind-difference method for aerodynamic problems governed by the Euler equations. In Lectures in applied mathematics, Part 2, Vol. 22. (pp. 327–336).