Discrete Solution of the Breakage Equation Using Markov Chains

8248

Ind. Eng. Chem. Res. 2010, 49, 8248–8257

Discrete Solution of the Breakage Equation Using Markov Chains Muammer Catak, Nursin Bas,* Kevin Cronin, John J. Fitzpatrick, and Edmond P. Byrne Department of Process & Chemical Engineering, UniVersity College, Cork, Ireland

Analytical solution of population balance equations (PBEs) may be impossible except for some simple cases. In the literature there are a number of methods to solve PBEs including discrete methods, Monte Carlo simulation, and method of moments. In this paper, the Markov chain is presented as a discrete solution for a population balance equation of a breakage process for determining the particle size distribution (PSD) over time. The transition matrix P, which is the key operator of a Markov chain, is built using breakage equations. Thereafter, from calculating transition matrix, P, the particle size distribution of the system is easily evaluated using the Markov chain. According to simulation results, if the size range of the system is divided into a sufficient number of states and an appropriate transition time step was chosen, then results from the Markov chain are in agreement with the analytical solution of PBEs governed by the same breakage functions. In addition to theoretical illustration, the Markov theory was employed to model the breakage process of aggregated food products passing through a pneumatic conveying pipeline rig. 1. Introduction Breakup process is an important topic in particulate systems including granulation, crystallization, grinding, and conveying. Particle breakage occurs due to particle-particle collision forces or particle-equipment interactions.1 Particle breakage can be a problem during conveyance, particularly if the particles involved are granular and/or friable. In general terms two breakage mechanisms for dry granules have been proposed; first erosion or attrition and second fracture or fragmentation.2 Where erosion is the dominant breakage mechanism there results one large fragment of size close to the parent aggregate and a number of smaller fine particles. Where breakage is by fracture, this results in the production of a number of smaller fragments. In addition, the fracture type of breakage is divided into two modes; cleavage, in which parent particles break into a small number of fragments of similar size and shattering which results in many fragments over a wide range of sizes.3 Particle breakage is usually considered an undesirable process since it can result in reduction in particle size, in changes to the particle size distribution, dust generation, and handling and storage problems which may cause the particles to not satisfy the requirement specifications any longer.4 On the other hand, it can be desirable and intentional depending on the application. For instance, particle is reduced intentionally during grinding and milling applications. Whether it is desirable or unwanted, mechanisms of the particle breakage should be investigated in order to control the process. Mathematical modeling tools are very helpful to tease out the background of the breakage. In literature, population balance modeling is commonly employed to model a breakage process.4,5 Population balance equations (PBEs) are underpinned by the law of mass conservation. A continuous density function in PBEs makes the model quite powerful in analyzing the dynamics of a process. However, the structure of these equations is complex due to the intrinsic partial integro-differential equations and hence analytical solutions may not often be possible except for rather simple cases. Nonetheless there are a number of numerical methods described in the literature which allow for PBEs to be solved.6–17 These methods include Monte Carlo * To whom correspondence should be addressed. Tel.: +353 214903096. E-mail: [email protected].

simulation, method of moments, successive approximations, and discrete formulations. Berthiaux18 has reported an application of the Markov chain on particle grinding processes. The author first introduced the Markov chain and concretized it using examples. Batch and continuous grinding processes are then analyzed by the Markov chain. The link and the consistency between the stochastic (Markovian) and deterministic (population balances) approaches are also discussed. Important in this paper is not just the linkage between the Markov theory and the population balances but also a demonstration of the possible cooperation of these two strong modeling approaches. The population balances approach is a powerful way to grab the basic mechanisms of the processes since its breakage functions are based on probabilistic theory. On the other hand, the Markov theory has the well-developed features of the population balances, while it is flexible such as Monte Carlo simulation, and the general behavior of the process can be captured at the same time with small details. Moreover, modeling of the particle breakage using the Markov chain is examined by Berthiaux et al.19 as well as the other main particulate processes. The multidimensional Markov chain modeling was explained in addition to the one-dimensional one. In this paper, a Markov chain model is developed to analyze the breakage process of aggregate food products during pneumatic conveying. In addition, it was shown that the Markov model might adequately be employed to the breakage equation as a discrete solution. 2. Markov Chains Consider a stochastic process {Xn; n ) 1, 2, 3, ...} that takes a countable or finite number of possible values, where Xn notates a random variable. If Xn ) i, then the process is said to be in state i at time n. It is supposed that whenever the process is in state i, there is a probability pij that it will next be in state j. It can be explained as20 pij ) P{Xn+1 ) j|Xn ) in, Xn ) in-1, ..., X1 ) i1} ) P{Xn+1 ) j|Xn ) in} (1) for all states i1, i2, ..., in and all n g 1. Such a stochastic process is known as a Markov chain. Equation 1 might be interpreted

10.1021/ie100216g  2010 American Chemical Society Published on Web 07/27/2010

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010

as stating that the conditional distribution of any future state of Xn+1, given past states X1, X2, ..., Xn-1 and the present state Xn, is independent of the past states and depends only on the present state. This is called as the Markovian property.21 Some constituents in certain observation processes (for instance particles) are experiencing a change in their properties (size, temperature, and location). These changes frequently have random characters, at least have some random components. For the aim of this work, the entire particle size range is divided into n finite intervals and a randomly selected particle belongs to one of these n intervals. In other words, the particle size range is converted to discrete values and can be represented with integers. The set of these discrete values of particle property is called state space. The probabilistic values of particle property over the state space constitutes the state probability Vector a with components ai(t) where i ) 1, 2, ..., n. a(t) ) [a1(t), a2(t), ..., an(t)]1×n By focusing on the single property of the particle, and then dividing the total process time into the finite periods of time τ, a representation of the process which is discrete in time as well as the particle property is obtained. The transition matrix P has entries pij which represents transition probabilities from state i to state j at each time step. On this basis, the transition matrix P for n states is given by an n × n matrix as follows:22

(

p11 p21 P) ... pn1

p12 p22 ... pn2

... ... ... ...

p1n p2n ... pnn

)

n×n

Properties of the transition matrix are (i) The sum of all probabilities of state i equals 1, that is, pi1 + pi2 + ... + pin ) 1, ∀ i ) 1, ..., n. (ii) Each row of the transition matrix has at least one nonzero element. (iii) Entities of the transition matrix are nonnegative, since pij is a probabilistic ratio, that is, pij values lie between 0 and 1. If the probability for an entity currently in the state j at time t is denoted by aj(t), then the state probability distribution of state i for the next time step, ai(t + τ), is given by the sum of product of all probabilities. This is formulated as follows: n

ai(t + τ) )

∑ p a (t)

(2)

ji j

j)1

Accordingly, eq 2 can be represented in matrix notation as

(

[a1(t + τ), ..., an(t + τ)] p11 p21 ) [a1(t), ..., an(t)] ... pn1

p12 p22 ... pn2

... ... ... ...

p1n p2n ... pnn

)

(3)

If the transition probabilities are independent of time, in other words the transition matrix P is a constant matrix, then the it is called a homogeneous Markov chain. Random walks is the most common example for the homogeneous chain.23 The game of Snakes and Ladders which has always been an eye catching example of Markov chains is another example for homogeneous chains.24 In chemical engineering, residence time distribution (RTD) can be modeled as a homogeneous Markov chain.25,26 On the other hand, if the transition matrix P is not constant over time, it can be either a linear nonhomogeneous or a nonlinear nonhomogeneous chain depending on the changing

8249

of P over time. The general form of eq 2 for time-homogeneous and nonhomogeneous transition matrices P for different time steps can be given by the following equations, respectively: a(t + mτ) ) a(t)pm

(4)

a(t + τ) ) a(t)p(t)

(5)

Although the Markov theory has been used in many fields after Andrei Andreyevich Markov released his work,21 it has rarely been applied in the chemical engineering area. There is only one book (which was published in 1998) on the application of Markov chains in chemical engineering.27 L.T. Fan and his colleagues were a source of Markov chains study from the 1970s to 1990s. Currently Henri Berthiaux is at the center of another group who work on Markov chains. A list of applications of Markov chains in chemical engineering and powder technology is as follows: particle flow in fluidized bed,28–30 particle mixing,31–35 particle breakage,18,19,36 deep-bed filtration,37 fluidized bed granulation.38 3. Population Balance Modeling of Particle Breakage Population balance modeling of a breakage process can be constructed using breakage functions. It might be either in continuous or in discrete form. In addition, two types of breakage equations are defined: (i) number-based breakage equation and (i) mass-based breakage equation. When particles display homogeneous density then a mass-based breakage equation can be considered as a volume-based breakage equation. 3.1. Continuous Population Balances. The general continuous-size form of the population balance equation for modeling a breakage process of a batch system can be given as follows:39 ∂f (x, t) ) ∂t

∫

∞

x

b(x, y) S(y) f(y, t) dy - S(x) f (x, t)

(6)

where x is particle size (equivalent diameter or volume) t is time, f (x,t) is the number density function of particles of size x at time t; b(x,y) is the probability distribution of particles of size x resulting from the breakup of particles of size y. S(x) is called the breakage frequency (or selection function) and is the rate at which particles of size x break per unit time. Equation 6 states that, the change in the number of particles of size x over an incremental time step at time t depends on the number of new particles of size x produced by break-up of a particle bigger than size x (the integral part in right-hand side) and depends on the average number of particles lost by breakage of particles of size x. Note that b(x,y) should satisfy the mass conservation requirement such that the total mass of particles resulting from the break-up of a particle of size of y must be equal to the mass of y.40 Since the procedure is same for both number-based and massbased breakage equations, only the number-based case will be detailed in the following sections and the equivalent particle diameter is taken to represent the corresponding size interval. 3.2. Discretized Population Balance. Classification of a continuous size range into discrete intervals is the first step in constructing a discretized population balance equation. Two main approaches can be used to select the classification of the states. (1) Uniform discretization: According to this approach, the interval of each state is constant. If the particle size is taken as its equivalent diameter, then r ) di - di-1.

8250


(2) Geometric discretization: In this approach, the interval of the next state is directly proportional to its previous state. If the particle size is its equivalent diameter, then di di-1

r)

where di-1 is lower limit of interval i, and di is upper limit of interval i. The underlying basis of discretization is to convert the main population balance of integro-differential equations into ordinary differential equations. Thus, the discrete-size form of a breakage equation can be given as dNi(t) ) dt

n

∑ b S N (t) - S N (t) ji j j

(7)

i i

{

∫

bji ) 0

di

di-1

b(x, dj-1) dx, if 2 e j e n and 1 e i < j

(8)

else

The common assumption in the literature is that if a particle of interval i breaks, its fragments cannot remain in the interval i, that is, the size of any fragment particle must be sufficiently large to produce a transition within the discretization scheme. The maximum interval of the fragment particle can be i - 1. Additionally, particles in the lowest interval are assumed to be the smallest possible size and they cannot break. On this basis, the breakage frequency of the first interval S1 equals 0.

Assuming that S(x) and b(x,y) are not affected by time, PBEs can be solved using a time-homogeneous Markov chain. It is also possible to apply different stationary (time-homogeneous) transition matrices for countable subtime intervals. The Markov chain algorithm can be applied to eq 7 by taking the Markovian time step as dt ) τ which yields the following equation: n

∑ b τS N(t) - τS N (t) ji

j

(9)

i i

j)i+1

In this equation, ∆Ni(t) is the rate of change in the number distribution of state i. When Ni(t) is added to both sides of the equality in eq 9, then the following equation is obtained: n

Ni(t) + ∆Ni(t) )

1 0 0 1 - τS2 . . D) . . 0 0 0 0

∑ b τS N (t) + (1 - τS )N (t) ji

j j

i

i

0 0 . . 0 0

. . . . . . . . . 1 - τSn-1 1 . .

(10)

j)i+1

The summation of the number distribution is represented by eq 10 and the rate of change in number distribution of interval i in a present state is equal to its next future state. That is to say that the left-hand side of the equality in eq 10 denotes the next future state of a given present state. 4.1. Transition Matrix P. The transition matrix P is a key element of the Markov chain model and construction of this matrix P for a breakage process is based on the characterization of a diagonal matrix D and a lower triangular matrix L. In this

0 0 . . 0 - τSn

)

(11)

n×n

2. Define an n × n lower triangular matrix L with entries: Li ) τSiQi

(12)

where Li represents ith row of L matrix, Qi is ith row of matrix Q which is built by qij (see eq 8), and Si is ith element of row vector S. In a breakage event when a particle breaks, each fragment of this particle can only go to lower states. Therefore the lower triangular matrix can be chosen since the upper entries are zero.

(

0 0 τS2b21 0 τS3b31 τS3b32 L) . . . . τSnbn1 τSnbn2

. . . . . .

. 0 . 0 . 0 . . . . . τSnbn,n-1

0 0 0 . . 0

)

n×n

According to eq 11, when selecting τ, it should be considered that τSi cannot be greater than 1. The Markovian transition matrix P equals sum of matrices L and D: P)L+D

4. Application of the Markov Chain to a Breakage Equation

∆Ni(t) )

(

Dii ) 1 - τSi

j)i+1

where Ni(t) is the number of particles in the interval i at time t, Si is the breakage frequency of particles in the interval i, and bji is the probability distribution of particles of interval i which are formed from a break-up of particles of interval j. Accordingly bji can be formulated as bji )

manner, eq 10 can be represented in matrix notation by dividing the whole system in two parts: 1. Define an n × n diagonal matrix D with entries:

(13)

Comments on Transition Matrix P. If number distribution is being modeled rather than mass distribution and the total number of particles increases in a breakage process, then the usual sum to unity will not apply directly. Initially the number distribution function, f (x,t), is normalized (i.e., area under f (x,0) is 1). However, the total number in the system will increase with time as a result of the breakage events. This results in the area under the curve, f (x,t), to be greater than 1 when t > 0. This is analogous to that for the transition matrix P. Even if the sum of initial state probability distribution vector, a(0), is 1, this sum will be greater than 1 at later stages. For the comparison, both f (x,t) and a(m) (m is a scaler which equals t/τ) should be normalized to 1 for each time step. On the other hand, there is no need for the renormalization procedure for the modeling of the mass distribution since the particle mass is an additive property. This is same for the transition matrix P which meets the condition of normalization under the case of mass distribution. 5. Results and Discussion 5.1. Case Study 1. In this section a realistic case study is performed for a binary breakage process, which means that a single particle gives two fragments, in order to compare results of the analytical solution which is obtained from eq 6 with results from the Markov chain model. 5.1.1. Analytical Solution. Let the specific breakage frequency be proportional to a characteristic diameter of the particle:


S(x) ) kx

3

8251

(14)

and define the fragment size distribution which is uniform with respect to volume size as11 bV(V,w) )

2 w

(15)

where k is a constant, w is the volume of parent particle, and V is the volume of the fragment. The physical meaning of eq 15 is that if a particle is broken the consequent fragments (daughter particles) can go into any smaller size intervals with an equal chance. To obtain a dimensionless fragment size distribution, dimension of the constant coefficient in eq 15 (that is 2) is chosen as the dimension of the parent particle’s volume. The relationship of the fragment size distribution between the particle volume and the equivalent diameter with the assumption of spherical particles can be defined as11 b(x,y) )

π 2 x b (V,w) 2 V

(16)

where V)

π 3 x 6 Figure 1. Fragment size probability distributions for different parent sizes.

and w)

π 3 y 6

As a result, the fragment size distribution with respect to the particle diameter is represented as b(x, y) )

6x2 y3

(17)

where y is the diameter of the parent particle and x is the diameter of the fragment. The fragment size probability distribution in eq 17 for the particle size range [0 mm, 2 mm] is shown in Figure 1. The PBEs for this breakage process with the breakage frequency in eq 14 and the fragment size probability distribution in eq 17 is written as ∂f (x,t) ∂t

∫

∞

x

6x2 3 ky f (y, t) dy - kx3f (x, t) 3 y

(18)

The analytical solution of eq 18 for the initial condition at t ) 3 0, f (x,0) ) 3x2e-x was found as41 f (x, t) ) 3x2(1 + kt)2e-x (1+kt) 3

(19)

The solution of the breakage process with the breakage frequency in eq 14 and the fragment size distribution in eq 17 will be carried out in the next section by applying the Markov chain to the same initial condition of eq 19. 5.1.2. Markov Chain Model. In this section an application of the Markov chain is introduced through the example that is considered in section 5.1.1. First, the particle size range should be divided into an appropriate number of states. The particle size range from 0 mm to 2 mm is divided into eight states by applying the uniform discretization that is given in section 3.2. The resultant values for these eight states and the corresponding initial particle size distribution are illustrated in Figure 2. Eight states is just an arbitrary number; it is chosen for convenience to illustrate the method succinctly. Second, the breakage frequency which is given in eq 14 should be described in discrete

Figure 2. Initial particle size distribution of the size range [0 mm, 2 mm] with representative average sizes of intervals.

form so that it is assigned for each state. The corresponding equation for the breakage frequency for a constant k can be defined as

{

Si ) kx3i if 2 e i e n Si ) 0 else

(20)

where xi is the representative diameter of state i. Thus, the discrete breakage frequency vector S can be written as S ) [0, S2, ..., Si, ..., Sn] The values of Si in eq 20 are calculated for each state of eight-state space by using the average particle sizes shown in Figure 2 and for a conjectural constant k ) 10-7 m-3 s-1 as

8252


S ) [0 0.005 0.024 0.067 0.142 0.260 0.429 0.659 ] (21) The diagonal matrix D, lower triangular matrix L, and transition matrix P are obtained by eqs 11-13 with the values in eq 21 and taking transition time step τ ) 1 s. The calculated results for matrices D, L, and P are shown as follows:

( ( (

D) 1 0 0 0 0 0 0 0 0 0.995 0 0 0 0 0 0 0 0 0.976 0 0 0 0 0 0 0 0 0.933 0 0 0 0 0 0 0 0 0.858 0 0 0 0 0 0 0 0 0.740 0 0 0 0 0 0 0 0 0.571 0 0 0 0 0 0 0 0 0.341

0 0 0.043 0.035 0.031 0.029 0.028 0.027

0 0 0 0.094 0.085 0.079 0.076 0.073

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.165 0 0 0 0.154 0.254 0 0 0.147 0.242 0.362 0 0.142 0.235 0.350 0.488

0 0 0 0 0 0 0 0

8×8

(23)

P) 1 0.105 0.006 0.005 0.005 0.004 0.004 0.004

8×8

(22)

L) 0 0.105 0.006 0.005 0.005 0.004 0.004 0.004

) ) )

0 0.995 0.043 0.035 0.031 0.029 0.028 0.027

0 0 0.976 0.094 0.085 0.079 0.076 0.073

0 0 0 0.933 0.165 0.154 0.147 0.142

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.858 0 0 0 0.254 0.740 0 0 0.242 0.362 0.571 0 0.235 0.350 0.488 0.341

8×8

(24)

The initial state probability distribution is found for the initial condition in eq 19 by the following equation: ai(0) )

∫

di

di-1

f (x,0) dx

(25)

Thus, the initial probability distribution for eight states by using the size ranges that are shown in Figure 2 is obtained as follows: a(0) ) [0.016 0.102 0.227 0.288 0.226 0.108 0.029 0.004 ] (26) 5.1.3. Discussion. The analytical solution obtained in eq 19 for time up to t ) 100 s and for k ) 10-7 m-3 s-1 is displayed in Figure 3. In this figure, the mean size of particles is monotonically decreasing with increasing time. The mean size is 0.9 mm at the start, but it decreases to 0.5 mm after 50 time steps as might be expected from a breakage process. In addition, the variance of the number distribution decreases with time. Figure 4 shows that the agreement between the analytical solution and the Markov chain with eight states is not satisfactory. Furthermore, it is clear that the percentage error between the analytical solution and the Markov chain increases with an increase in the number of time steps (see Table 1). The

Figure 3. 3D visualization of analytical solution for size range [0 mm, 2 mm] and time up to 100 s.

percentage error is obtained from the ratio of the absolute difference between the analytical solution and the model to the analytical solution multiplied by 100. After the fifth time step there is a 30% error. The two most probable reasons for the difference are (i) the number of states, i.e., 8, is not enough to represent the system; thus the number of states should be increased to get more accurate results. (ii) As introduced in section 3.2, discretization may produce some intrinsic problems. For example, suppose there are particles ranging from 0 to 90 unit volume and they are classified into five intervals as [0 10], [10 - 30], [30 - 50], [50 - 70], [70 - 90] unit volume. By considering only binary breakage, if a particle has 80 unit volume and fragments into sizes of 5 and 75 unit volume, then there is no change in the number of particles in state 5 ([70 90] unit volume). To overcome this problem, discretization of size domain should be well-defined. A large number of states reduces the likelihood of these difficulties. By using software, such as Matlab, the Markov chain with several thousands of states can be easily implemented. It could be expected that the use of 100 states will give better results than 8 states. On this basis, a 100-state space model has been simulated. The state space is obtained by dividing the particle size range [0 mm, 2 mm] into 100 states via the uniform discretization similar to eight-state space. The Markov chain simulation is compared with the analytical solution in eq 19 for the certain times. Figure 5 displays the simulated results of the Markov chain and analytical solution for the 100-state space. It can be seen that for the same transition time (τ ) 1 s); increasing the number of states results in a more accurate solution. For the eight-state space (see Figure 4), the agreement between both results, the Markov chain and the analytical solution, is satisfactory up to t ) 10 s. However, this agreement becomes deficient with time. The overall accordance between the Markov chain and the analytical solution are comparatively sufficient for the100-state space model shown in Figure 5. Consequently, it is obvious that increasing the number of states in the system leads to more accurate results. To support this idea, the Markov chain has been simulated for 20- and 60state spaces as well. The simulation results, especially for the 100-state space model, shows that the Markov chain gives a very good agreement with the analytical solution, as illustrated in Figure 5. When the whole size domain is divided into enough number of states, then the Markov chain can represent the analytical solution with high accuracy.


8253

Figure 4. Simulations of the Markov chain model and the analytical solution for eight-state space model. The agreement between both results becomes less with time. Table 1. Percentage Errors between the Analytical Solution and the Markov Chain Models According to Number of State Space percentage error (%) time [s]

8 states

20 states

60 states

100 states

1 5 10 50 100

13.1 30 39.9 45.7 48.8

3.5 7.1 10.2 14.3 19.8

0.4 1.1 1.9 2.6 3.3

0.2 0.8 1.2 2.3 2.7

Furthermore, statistical analysis of the results was also undertaken. The mean sizes for various time steps are calculated using the following equation: n

xmean(t) )

∑ x a (t) i i

(27)

i)1

Similarly, the variance in particle size is obtained from n

Var(t) )

∑ a (t)(x i

i

- xmean(t))2

(28)

i)1

The standard deviation denoted by σ, is the square root of the variance.

Comparisons of results obtained by the Markov chain and the analytical solution for the mean of the number distribution at different time steps are displayed in Figure 6. Mean size prediction by the Markov chain converges to the analytical solution with an increase in the number of state space. On the other hand, the standard deviation in particle size prediction of the Markov chain is very similar to the analytical solution for the different number of state space as illustrated in Figure 7. A probability distribution function is characterized by shape, location, and scale parameters. For the normal distribution, for instance, location and scale parameters correspond to the mean and standard deviation, respectively. For a distribution, the location parameter is dominant on the mean, while the scale and the shape parameters mostly affect the standard deviation. According to Figures 6 and 7, the Markov chain has good agreement with the analytical solution with respect to the scale and shape parameters even for a small number of states. However, for the location parameter accuracy increases with the number of state space. Fractional error, which is found by the ratio of the absolute difference between the analytical solution and the model to the analytical solution, of the models versus number of state space at t ) 100 s is shown in Figure 8, where the larger number of states produces more accurate

8254


Figure 5. Simulations of the Markov chain model and the analytical solution for 100-state space model. The agreement between both results is efficient with time.

Figure 6. Comparison of the mean sizes for the analytical solution and for various state space models of the Markov chain model with time.

Figure 7. Comparison of the standard deviation over size for the analytical solution and for various state space models of the Markov chain model with time.

results. In contrast to the mean, the standard deviation is not greatly influenced by number of state space. 5.2. Case Study 2. In this section, an application of Markov chain to a real breakage process is presented. The related experimental data is taken from literature.42 They applied the population balance model to characterize the breakage process of aggregated food products passing through a pneumatic conveying pipeline rig. In the experiment, the granola products were formed subject to impeller agitation at 300 rpm for 9 min using a high shear mixer, then the aggregates were baked in an oven for 10 min at 170°. Trials were carried out by applying compressed air driving pressures of 200, 300, and 400 kPa which

resulted in flow velocities of 23, 34, and 42 m s-1, respectively. The particles were recycled through the rig a number of times, 1, 2, 3, 4, 5, 10, 15, 20, and the particle size distribution of the samples were recorded. In their study, they assumed that particles are spherical and of broadly homogeneous density, thus a mass-density-based breakage equation was used to model volume fraction. Particle breakage mainly occurs during conveyance as product is transferred on its way to packaging due to the particle-wall collisions. In this section, breakage of the same process is modeled using the Markov chain, and the results obtained from the Markov model and the experimental results are compared.


Figure 8. Fractional error of the simulation regarding to number of state space at t ) 100 s.

5.2.1. Markov Chain Model. The experimental data attained by Bas et al.42 range from 0.5 mm to 11.62 mm with the number of intervals equal to 30. The particle diameter was discretized by using the uniform discretization with a constant difference r ) 0.4 mm. Additionally, a representative size of each interval is calculated using arithmetic mean of size ends. Thus, a priori 30-state space is used for the Markov chain. Moreover, experimentally given initial size distribution is taken as the initial state probability distribution of the Markov chain. It is assumed that the breakage frequency per cycle is proportional to its representative diameter, then discrete breakage frequency is given in the following form: Si ) kxi

3

(29)

where xi is the representative particle diameter and k is the breakage constant which is depended on the average air velocity in the pipe corresponding applied pressures. Berthiaux et al.19 used the following formulas for particle breakage in which the particle can transit only into neighboring intervals as bjm )

xn3R xj-13R + xm3R

xj3R x13R + x23R + ... + xj3R + ... + xi-13R

)

xj3R i-1

∑x

3R

k

k)1

(31) The value of R is selected as an adjustment parameter of the model, and it should be in the range 1 < R < 2 in applications to model a real process.19 The discretized volume-based breakage function bji is the volume fraction of particles broken from interval j that go to interval i. Any broken particle in interval j cannot go into interval j, that is, j-1

∑b

ji

i)1

)1

The transition matrix was built as explained in section 4 using the knowledge of discrete breakage functions. The model constants R and k were initially obtained by minimizing the error between predicted and experimental particle size distributions. For the validation of the Markov model, an independent set of experimental data obtained by following the same experimental procedure was used to refine these parameters. Using the maximum likelihood approach, which aims to achieve the minimum squared error between the model and the experimental measurements of the size distribution, the optimum R and k are taken as 1.07 and 7.5 × 10-11 m-3 s-1. The transition time step is taken as the corresponding number of cycles of the granola under different pressures. Consequently, all necessary elements such as transition matrix, initial state probability distribution, and transition time step are obtained to run the Markov chain simulation. 5.2.2. Discussion. Comparison of the experimental results and the Markov model at a pressure 400 kPa after various number of cycles are displayed in Figure 9. It is clear that the Markov chain has sufficiently captured the particle size evolution of the breakage process of aggregate food products over time. The Kolmogorov-Smirnov test was applied as a statistical analysis, and it was found that the experimental and the model results have a consistency not less than 95% for each cycling time. In addition, the mean particle size evolution, which is given in Table 2, is examined. The mean size, which is 5.18 mm at the initial, reduces in increasing number of cycles as expected from a breakage process. According to the results, the Markov chain is an efficient tool to predict the mean particle size as well as the overall particle size distribution. 6. Conclusion In this paper, the Markov chain theory is applied to the particle breakage process. First it was shown that if the size range of the system is divided into a sufficient number of states (n) and an appropriate transition time step (τ) is chosen, then the Markov chain exhibits a very good approximation to the analytical solution of corresponding population balance equation. In theory, if n approaches infinity while τ goes to zero, the Markov chain becomes identical to the analytical solution: lim lim {Markov chain} f analytical solution

(30)

If the restriction on only permitting passage to adjacent states is removed, then the distribution function in eq 30 can be generalized as bij )

8255

(32)

nf∞ τf0

However, in practice it is impossible. Therefore, a well chosen (n,τ) pair is necessary for an efficient Markov chain model. The sufficiency of the number of states and transition time step are dependent on the process kinetics.The authors suggest to divide the state property under analysis into 30 intervals as an initial trial. Some process systems have a natural periodicity such as rotating mixers (period of rotation) or vibratory conveyors (period of excitation) for which the transition time step can be chosen from basic physics. For systems which have no period, the time step should be chosen taking process kinetics and total processing time into account. In addition to the theoretical case study, an application of the Markov chain to a real breakage process is presented. The results show that the Markov model is an adequate tool to predict the particle size distribution over time. In literature, it has been shown that the Markov chain is a proper modeling way for grinding processes. In addition to this, the Markov chain might be employed as a useful tool for the solution of PBEs of the breakage process. The Markov model developed in this study is not a specific model for a certain breakage process. It is an advantage of the

8256


Figure 9. Volume fractions after a various number of cycles at a pressure of 400 kPa. Table 2. Predicted and Experimental Mean Sizes at a Number of Cycles for the Pressure 400 kPa and Corresponding Percentage Errors cycle

mean size (experimental) [mm]

mean size (model) [mm]

percentage error [%]

0 1 2 3 4 5 10 15 20

5.18 5.08 4.98 4.90 4.82 4.75 4.43 4.19 4.01

5.18 4.75 4.64 4.54 4.46 4.35 4.20 4.12 4.04

0 6.43 7.01 7.45 7.53 8.34 5.41 1.87 1.02

model to be applicable to any systems. On the other hand, it cannot give information of the processes in detail. For future work, nonhomogeneous (time-dependent) and/or nonlinear (state-dependent) Markov modeling of grinding processes might be developed. Such a Markov model may succeed where PBE failed in describing what is usually called “nonlinear grinding” in literature. Acknowledgment This publication has emanated from research conducted with the financial support of Science Foundation, Ireland. Authors

N. Bas, E. Byrne, and J.J. Fitzpatrick thank the National Development Plan, through the Food Institutional Research Measure, administered by the Department of Agriculture, Fisheries and Food, Ireland, for providing funding for this work. Literature Cited (1) Salman, A.; Ghadiri, M.; Hounslow, M. Handbook of Powder Technology: Particle Breakage; Elsevier; Amsterdam, London, 2007. (2) Iveson, S.; Litster, J.; Hapgood, K.; Ennis, B. Nucleation, growth and breakage phenomena in agitated wet granulation processes: A review. Powder Technol. 2001, 117, 3–39. (3) Redner, S. Fragmentation Statistical Models for the Fracture of Disordered Media; Elsevier: Amsterdam, 1990. (4) Salman, A.; Reynolds, G.; Hounslow, M. Particle impact breakage in particulate processing. KONA 2003, 21, 88–99. (5) Hounslow, M.; Pearson, J. M. K.; Instone, T. Tracer studies of highshear granulation: II Population balance modeling. AIChE J. 2001, 47, 1984– 1999. (6) Ramkrishna, D. Solution of population balance equations by the method of weighted residuals. Chem. Eng. Sci. 1971, 26, 1134–1136. (7) Austin, L. G. Introduction to Mathematical Description of Grinding as a Rate process. Powder Technol. 1972, 5, 1–17. (8) Bapat, P.; Tavlarides, L.; Smith, G. Monte Carlo simulation of mass transfer in liquid-liquid dispersions. Chem. Eng. Sci. 1983, 38, 2003–2013. (9) Ziff, R. New solutions to the fragmentation equation. J. Phys. A 1991, 24, 2821–2828. (10) Hounslow, M.; Ryall, R.; Marshall, V. A discretized population balance for nucleation, growth, and aggregation. AIChE J. 1988, 34, 1821– 1832.

Ind. Eng. Chem. Res., Vol. 49, No. 17, 2010 (11) Hill, P. J.; Ng, M. K. New discretization procedure for the breakage equation. AIChE J. 1995, 41, 1204–1216. (12) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-I: A fixed pivot technique. Chem. Eng. Sci. 1996a, 51, 1311–1332. (13) Kumar, S.; Ramkrishna, D. On the solution of population balance equations by discretization-II: A moving pivot technique. Chem. Eng. Sci. 1996b, 51, 1333–1342. (14) Vanni, M. Discretization procedure for the breakage equation. AIChE J. 1999, 45, 916–919. (15) Diemer, R.; Olson, J. A moment methodology for coagulation and breakage problems: Part 2-Moment models and distribution reconstruction. Chem. Eng. Sci. 2002, 57, 2211–2228. (16) Kostoglou, M.; Karabelas, A. An assessment of low-order methods for solving the breakage equation. Powder Technol. 2002, 127, 116–127. (17) Marchisio, D. L.; Pikturna, J. T.; Fox, R. O.; Vigil, D. R.; Barresi, A. A. Quadrature method of moments for population balance equations. AIChE J. 2003, 49, 1266–1276. (18) Berthiaux, H. Analysis of grinding processes by Markov chains. Chem. Eng. Sci. 2000, 55, 4117–4127. (19) Berthiaux, H.; Mizonov, V.; Zhukov, V. Application of the theory of Markov chains to model different processes in particle technology. Powder Technol. 2005, 157, 128–137. (20) Ross, S. M. Stochastic Processes; John Wiley and Sons: New York, 1996. (21) Markov, A. Extension of the law of large numbers to dependent events. Bull. Soc. Phys. Math., Kazan, Russ. 1906, 15, 155–156. (22) Farina, L.; Rinaldi, S. PositiVe Linear Systems: Theory and Application; Wiley: New York, 2000. (23) Kemeny, J. G.; Snell, J. L. Finite MarkoV Chains; D. Van Nostrand: New York, 1960. (24) Cronin, K.; Bas, N.; Byrne, E. Development of a computer-based teaching tool in probabilistic modelling. International Symposium for Engineering Education, Dublin City University, Ireland, 2007. (25) Fan, L. T.; Too, J. R.; Nassar, R. Stochastic simulation of residence time distribution curves. Chem. Eng. Sci. 1985, 40, 1743–1749. (26) Berthiaux, H.; Mizonov, V. Applications of Markov chains in particulate process engineering: A review. Can. J. Chem. Eng. 2004, 82, 1143–1168. (27) Tamir, A. Applications of MarkoV Chains in Chemical Engineering; Elsevier: Amsterdam, The Netherlands, 1998. (28) Fox, R.; Fan, L. Stochastic analysis of axial solids mixing in a fluidized bed. Chem. Eng. Commun. 1987, 60, 27–45.

8257

(29) Dehling, H.; Hoffmann, A.; Stuut, H. Stochastic models for transport in a fluidized bed. SIAM J. Appl. Math. 1999, 60, 337–358. (30) Harris, A.; Thorpe, R.; Davidson, J. Stochastic modelling of the particle residence time distribution in circulating fluidised bed risers. Chem. Eng. Sci. 2002, 57, 4779–4796. (31) Wang, R. H.; Fan, L. T. Stochastic modeling of segregation in a motionless mixer. Chem. Eng. Sci. 1977, 32, 695–701. (32) Fan, L. T.; Shin, S. H. Stochastic diffusion model of non-ideal mixing in a horizontal drum mixer. Chem. Eng. Sci. 1979, 34, 811–820. (33) Aoun-Habbache, M.; Aoun, M.; Berthiaux, H.; Mizonov, V. An experimental method and a Markov chain model to describe axial and radial mixing in a hoop mixer. Powder Technol. 2002, 128, 159–167. (34) Berthiaux, H.; Marikh, K.; Mizonov, V.; Ponomarev, D.; Barantzeva, E. Modeling continuous powder mixing by means of the theory of Markov chains. Part. Sci. Technol. 2004, 22, 379–389. (35) Ponomarev, D.; Mizonov, V.; Gatumel, C.; Berthiaux, H.; Barantzeva, E. Markov chain modelling and experimental investigation of powdermixing kinetics in static revolving mixers. Chem. Eng. Process. 2009, 48, 828–836. (36) Chakraborty, J.; Ramkrishna, D. Identification of Markov matrices of milling models. Ind. Eng. Chem. Res. 2009, 48, 9763–9771. (37) Nassar, R.; Chou, S. T.; Fan, L. Modelling and simulation of deepbed filtration: A stochastic compartmental model. Chem. Eng. Sci. 1986, 41, 2017–2027. (38) Catak, M.; Bas, N.; Cronin, K.; Tellez-Medina; D.; Byrne, E. P.; Fitzpatrick, J. J. Chem. Eng. J., in press, doi:10.1016/j.cej.2010.02.022. (39) Ramkrishna, D. Population Balance: Theory and Applications to Particulate Systems in En-gineering; Academic Press: San Diego, CA, 2000. (40) McGrady, E.; Ziff, R. M. “Shattering” transition in fragmentation. Phys. ReV. Lett. 1986, 58, 892–895. (41) Patil, D. P.; Andrews, J. R. G. An analytical solution to continuous population balance model describing floc coalescence and breakagesA special case. Chem. Eng. Sci. 1998, 53, 599–601. (42) Bas, N.; Pathare, P.; Catak, M.; Fitzpatrick, J.; Cronin, K.; Byrne, E. Mathematical modelling of granola breakage during pipe pneumatic conveying. Powder Technol., in press, doi:10.1016/j.powtec.2010.06.015.

ReceiVed for reView January 29, 2010 ReVised manuscript receiVed July 4, 2010 Accepted July 8, 2010 IE100216G

Discrete Solution of the Breakage Equation Using Markov Chains

Recommend Documents