Fast Compartmental Monte Carlo Simulation of Population Balance

Article pubs.acs.org/IECR

Fast Compartmental Monte Carlo Simulation of Population Balance Models: Application to Nanoparticle Formation in Nonhomogeneous Conditions Roberto Irizarry* Electronic Technologies, DuPont, 14 T.W. Alexander Drive, Research Triangle Park, North Carolina 27709, United States ABSTRACT: A new compartmental Monte Carlo (CMC) algorithm is introduced for the stochastic simulation of population balance models in spatially heterogeneous systems. The heterogeneities are modeled using a network of compartments. The algorithm is based on a new stochastic procedure called particle bundle flow (PBF) to model the stochastic transfer of particles between compartments in a given time interval (a time-driven algorithm). Different from other time-driven methods, the accuracy of the PBF is independent of the particle concentration. The validity of the PBF method is demonstrated by construction and confirmed with numerical experiments. A new strategy for time step control is developed to set bounds on the calculation of the time steps during the simulation. The CMC algorithm, based on the combination of the PBF algorithm with the τ point ensemble Monte Carlo algorithm, is a general-purpose methodology that can be applied to any network of compartments. The computational speed and the low computational load of this algorithm allow the solution of problems that may be intractable otherwise. A new hybrid strategy for the solution of problems with stochastic fluctuations and disparate time scales is also developed in this work. The CMC is applied to study the formation of nanoparticles in a large reactor utilizing a two-compartment model. The Monte Carlo ability to track single events is utilized to study the impact of turbulence and the stability factor on the generation of large particles. fluid flow cannot be ignored. The development of such methods is a very significant tool in the analysis of the scale-up and the development of optimal industrial operating conditions. Only recently have attempts been made to couple MC simulations with computational fluid dynamics (CFD).7 This type of rigorous approach is very effective, but for modeling the full volume of an industrial reactor can be computationally intensive. For example, in industrial reactors the number of CFD cells and number of simulation particles may become very large. A typical industrial reactor for the precipitation of metal powders will have on the order of 1021 particles, while a typical homogeneous MC simulation would have 103 to 105 simulation particles. Even for the case when direct numerical solution of a PBM can be coupled with a CFD model, computationally taxing simulations result.8 To reduce the complexity and computational cost of detailed heterogeneous simulation, one approach is to model the system with a low-dimensionality network of well-mixed compartments. Although a compartmental model is a very simple and traditional type of model, there are many new methodologies to break down complex systems into such simple compartmental models. These types of models have proven to be effective in reproducing the main dynamics of the detailed model with a reduced dimensionality and computational cost. Several heuristic rules have been presented in the literature to set up the compartment model network from CFD simulations.9−12 Knowledge of the scaling parameters and flow structure from finite element simulation has been used to develop compartment models as well.13,14 In addition to capturing nonhomogeneities in reactors, compartmental modeling has many applications in many heterogeneous complex

1. INTRODUCTION Stochastic simulations are important tools in the modeling of complex phenomena. One example of these types of simulations is the Monte Carlo (MC) solution of complex population balance models (PBMs).1,2 These are particle-based methods where the history of a set of particles in a simulation volume is followed as a function of time. Acceleration techniques for MC simulation of population balance models have been developed recently on the basis of a coarse-grained Markov model that approximates the actual PBM model.3−6 Although all these MC methods are very powerful, an important constraint of these methods is the assumption of a well-mixed system. In many situations the system cannot be approximated as spatially homogeneous. In these cases the PBM is coupled with the fluid flow equations. One example of this type of problem is the particle formation in industrial size reactors, as shown schematically in Figure 1a.

Figure 1. Particle formation coupled with fluid flow: (a) detailed flow field, (b) compartmental approximation of the flow field.

Received: Revised: Accepted: Published:

Less attention has been paid to developing MC simulation strategies to solve these types of problems, where the effect of © 2012 American Chemical Society

15484

April 28, 2012 August 28, 2012 October 23, 2012 October 24, 2012 dx.doi.org/10.1021/ie3011116 | Ind. Eng. Chem. Res. 2012, 51, 15484−15496

Industrial & Engineering Chemistry Research

Article

systems, including biological systems,15,16 system-based simulations, and particulate processes,17−19 among other applications. In general, there are two types of compartmental models. One type is shown in Figure 1b, where the flow field is replaced with a compartment network that has no direct connection with an actual physical space. The second type of model is shown in Figure 2, where in this case the actual domain is divided into

Figure 2. Particle formation coupled with fluid flow: (a) coarse partition of the volume, (b) compartments from the coarse partition. Figure 3. Schematic of the time-driven compartmental Monte Carlo algorithm.

coarse elements (Figure 2a) that become the compartments in the simplified model (Figure 2b). There are two categories of MC simulations, time-driven and event-driven simulations. The event-driven simulations are exact methods where a trajectory is generated by sampling the time for the next event (i.e., when the next collision will occur) and the actual event that will occur (i.e., which two particles will collide in this event). This process is continued until a final time is reached. In this exact simulation, the time and the events are stochastic variables. Although these are exact methods, this one event at a time approach can be very slow in many instances. To speed up the simulation, the time-driven simulation is based on setting a larger time step and then firing many events in that time interval. In this type of simulation, the time interval is a deterministic variable (fixed by the user) and the number of events and the actual events are the stochastic variables being sampled from probability distributions. The time-driven simulations are much faster than the event-driven simulations; on the other hand, these are approximated methods valid when certain conditions are satisfied (i.e., high concentration). When these conditions are not satisfied, the event-driven algorithm is the only alternative. One example of this situation is the simulation of very dilute systems, where the particle properties are not a deterministic number but a fluctuating variable with a probability distribution function (PDF). In this case, the exact event-driven simulations can generate the correct PDF by running a large number of simulations to generate the statistics. This work has three objectives. The main objective is to develop a fast and accurate general-purpose algorithm for the time-driven MC simulations of compartmental PBMs. The second objective is to develop a hybrid event-driven/time-driven algorithm to speed up the simulation of problems with fluctuating output that also have disparate time scales. The third objective is to apply the time-driven algorithm to simulate the formation of nanoparticles under turbulent reactor mixing. Because the algorithm is composed of several building blocks, a conceptual picture of the algorithm and its key features are briefly described in the following subsections. 1.1. Algorithm Summary. Figure 3 shows a conceptual schematic of the time-driven compartmental Monte Carlo algorithm (CMC). The algorithm starts with the assumption that the fluid flow is represented with a compartmental model. Another preprocessing

step is to set the upper bound for the maximum allowed time step. This bound will guarantee an accurate modeling of the particle transfer between compartments. The calculation procedure to set this upper bound is introduced in this work. After these prepossessing steps, the CMC algorithm integrates the model dynamics as follows: First, the time step for the current iteration is calculated. Then the time-driven MC step in each compartment is executed. Since each compartment is homogeneous, an existing time-driven MC method is used in this step. The next step is to select the particles to be transferred between compartments. This step is accomplished using a new procedure called the particle bundle flow (PBF) method. These steps are repeated until a final time is reached. The main objective of this work is the simulation in algorithm 3. A second algorithm designed to solve multiscale problems (to be discussed later) is also introduced in this work. 1.2. Algorithm Key Features and Relevance. The CMC proposed in this work is a new methodology to solve PBMs that are coupled with complex mixing and fluid flow conditions. The main advantage of this method comes from its low CPU time. This algorithm can be utilized to solve problems that were previously unfeasible or feasible only on high-performance clusters. This is a general-purpose methodology that applies to any type of compartmental model. The only requirement for the application of the algorithm is to specify the flows between compartments and that each compartment is treated as a well-mixed compartment. The PBF method is a new time-driven procedure developed in this work to model the transfer of particles from one compartment to another (in a network of interconnected compartments). This method is differentiated from other existing timedriven algorithms in the fact that its accuracy is independent of the number of particles in the compartment. This unique property of the PBF method can be utilized to accurately model very dilute compartments, including capturing the correct flow fluctuation statistics. This property enables the use of the PBF method to develop a novel hybrid strategy to solve problems with multiple time scales with fluctuating outputs in an efficient way. This hybrid strategy is a new algorithm introduced in this work that uses the PBF to accelerate the event-driven MC simulation. 15485

dx.doi.org/10.1021/ie3011116 | Ind. Eng. Chem. Res. 2012, 51, 15484−15496


Article

Different from the time-driven simulation, this hybrid algorithm can capture the statistics of small subpopulations. Other existing time-driven algorithms that use Poisson distributions (τ-leap methods) are based on the assumption of a large number of particles, limiting their application to systems with relatively high concentration (see Appendix C and references therein). Therefore, they cannot be used for the hybrid implementation. Furthermore, these methods will provide inaccurate simulation in dilute compartments. A new time step calculation procedure is introduced. This procedure sets an upper bound on the time step used in each iteration. This bound guarantees an accurate modeling of the particle flow. One advantage of this algorithm is that, once coded, it can be used for any network of compartments in a seamless way. The time-driven MC method used for the simulation inside each compartment is the τ point ensemble Monte Carlo (τ-PEMC) procedure.4 In principle, any existing MC method could be utilized in this step. The τ-PEMC method was selected for this task because it is a time-driven algorithm that can easily be interfaced with the PBF method. Furthermore, this algorithm accelerates the simulation speed by orders of magnitude as discussed in detail in refs 3 and 4. Thus, the combination of the PBF with the τ-PEMC method results in a very fast simulation methodology. This paper is organized as follows. The PBF utilized to model the particle transfer between compartments is introduced in section 2. The algorithm to set the upper bound on the maximum allowed time step is discussed in section 2.1. The event-driven PEMC and time-driven τ-PEMC methods are reviewed in section 3. The CMC simulation is discussed in section 4. The extension of this algorithm to solve muti-time-scale problems is discussed in section 4.1. The performance of the PBF method is discussed in section 5. The formation of metallic nanoparticles using the timedriven CMC is discussed in section 6. The conclusions are discussed in section 7.

PBD(K ; p , k max ) =

K max p K (1 − p)K max − K K! (K max − K )!

(1)

This binomial process only determines the number of particles leaving the compartment (Kout i ), but not the actual particles and their new compartment location. This is determined in a second step. Each of the Kout i particles is selected randomly with equal probability using a uniform distribution. For each selected particle, the new compartment j is selected with probability proportional to the flow intensity, f ij/f i. The next compartment j is therefore found by solving the following equation, such that j is the smallest integer satisfying it: rf i ≥ ∑jk = 1f ik, where r is a random number with uniform distribution in the interval [0, 1]. The algorithm is summarized next. Algorithm 1: particle bundle f low method, PBF(Δt) 1. For each compartment i: Sample Kout from a binomial i distribution, PBD(Kout ;p ,N ). i i c(i) 2. For each compartment i: Transfer Kout i random particles out of compartment i to a collection compartment (Ncomp + 1). For each particle: a. Select a random carrier particle, n, from compartment i. b. Find the next compartment j, where j is the minimum number such that rf i ≥ ∑jk = 1f ik, r being a random number from a uniform distribution. Remove particle n from compartment i and add it to a temporary accumulation compartment (Ncomp + 1). (Add the label of the found compartment, j, to the particle for future redistribution from compartment Ncomp + 1 to j: jnext(n) = j.) 3. Redistribute all particles n accumulated in Ncomp + 1 to their new compartment using the pointer jnext(n). This algorithm transfers temporarily all particles to a temporary compartment (Ncomp + 1) instead of to their new compartment, jnext(n). Sending the particles to this temporary compartment as an intermediary step eliminates the effect of any possible bias due to the compartment execution order. This is valid because a particle cannot hop more than one compartment (hopping constraint). In step 3, these particles accumulated temporarily in Ncomp + 1 are then sent to their final destination using the pointer jnext(n). This algorithm is faster and more accurate than other algorithms such as the time-discretized Markov chain method to simulate particle transfer.20,21 In this method, each particle has to be sampled, making the system computationally expensive if the system has a large number of particles. Also, the transition probabilities are determined by forcing the random time for the next event to be a deterministic fixed number (the selected fixed time interval). Binomial or Poisson processes have been used to model reaction−diffusion processes. In these models, the diffusion is considered as a first-order reaction or a reversible set of firstorder reactions22,23 (also see Appendix C). The method introduced in this work also results in a binomial process, but it is constructed by exploiting the statistical properties of well-mixed compartments, independent of the particle concentration (see Appendix A), instead of assuming first-order chemical reactions. This results in a process that is valid for compartments with very large or very small numbers of particles (a major departure from other τ-leap methods). In τ-leap methods, the time step has to be such that the number of species does not change much to satisfy

2. TRANSFER OF PARTICLES BETWEEN COMPARTMENTS: PBF METHOD A critical building block in the development of the CMC algorithm is an efficient and accurate procedure to simulate the transfer of particles between compartments. The method developed in this work is based on a construction of a binomial process for the number of particles leaving a well-mixed compartment. These steps are discussed in detail in Appendix A. By construction, this method is stochastically exact if the time step is selected such that the probability of a particle jumping more than one compartment in the given time interval is negligible. I call this constraint the hopping constraint. This is a time-driven method to transfer a bundle of particles in a given time interval. The method described next is valid for any concentration of particles (from very dilute to very concentrated). Consider a network of interconnected well-mixed compartments numbered 1, 2, ..., Ncomp with volumes Vi. The network description is completed with the specification of the steady flow rates between compartments, qij. Let f ij = qij/Vi and f i = ∑j f ij. Given a system with Nc(i) particles in compartment i, the objective is to determine how many particles will move from compartment i to a connected compartment j in a given time interval, Δt. If this time interval satisfies the hopping constraint, as shown in Appendix A, the number of particles leaving compartment i, Kout i , is a random number taken from a binomial distribution with parameters pi = f iΔt and maximum number of outcomes Nc(i): PBD(Kout i ;pi,Nc(i)), where 15486



Article

the τ constraint, that is, the propensities do not change during the time interval.24 2.1. Determination of Upper Bound UBh To Satisfy the Hopping Constraint. One intrinsic assumption in the PBF method is that the time interval is small enough that the particle cannot visit more than one compartment. This can be achieved by selecting a small enough time step. To maximize the time step utilized, the following algorithm can be used to find the upper bound UBh on the time interval that satisfies the particle hopping constraint to a given tolerance. The method consists of simulating a swarm of “test” particles distributed randomly in the network of compartments and following their history in a given time interval, Δt. The ratio of particles that visited more than one compartment in this interval to the total number of initial particles is utilized to estimate the probability of a particle visiting more than one compartment, P≥2(Δt). The procedure to determine P≥2(Δt) is described in algorithm 2. Enough test particles are used to guarantee that the calculation is accurate. Given a desired level of tolerance, εh (%) (typically 0.1−5%), UBh can be easily estimated by plotting P≥2(Δt) × 100 versus Δt and then locating the desired εh (%). Algorithm 2: determination of the percentage of particle hopping constraint violations 1. Input: Δt. 2. Set Nfail = 0. 3. Loop over Ntest test particles. a. Set Nvisits = 0, t = 0. b. Start the test particle in a random compartment i. c. Do while t < Δt: i. Determine the exit time, τi, by sampling τi = −ln(r)/f i, where r is a uniformly distributed random number. ii. Set t = t + τi. iii. Nvisits = Nvisits + 1. iv. Find the next compartment j, such that j is the smallest integer satisfying Rf i ≥ ∑jk = 1f ik, where R is a random number. v. Set j as the new current compartment (set i = j). End do while. d. If Nvisits > 2, then Nfail = Nfail + 1. End of loop over all test particles. 4. Output: P≥2(Δt) (%) = 100(Nfail/Ntest). In addition to the tolerance on the hopping constraints, another constraint that may be added is the percentage of transfers accepted in one time interval: Δt ≤ (α /max{fij }) ≡ UBf

on the simulation of the PERP model, which is a coarse-grained Markov model that approximates the original Markov model (the simulation box). The PERP model is described in section 3.1 for the case of nucleation aggregation, and the algorithms are described in sections 3.2 and 3.3. 3.1. PERP Model. First, the volume space is partitioned with M grid points, vi. All particles with volume close to a given grid point are collected in subensembles Φi called point ensembles. Let Ni be the number of particles in point ensemble Φi. The state vector is then defined as (N, Φ), where N = (N1, N1, ..., NM)T and Φ = (Φ1, Φ2, ..., ΦM)T. Here, each grid point is viewed as a chemical pseudospecies, Si, with Ni molecules in the simulation volume. For example, Figure 4 shows a partition of 17 particles (each with a set of different properties) according to their sizes

Figure 4. Schematic of the PERP jump Markov model.

into 5 point ensembles. The five pseudospecies (S1, S2, ..., S5) defined by this partition have (3, 4, 6, 3, 1) molecules in the simulation volume. In this model, each aggregation event, Eagg(i, j), between particles in Φi and particles in Φj consists of the following steps: Select random particles n and m from point ensembles i and j (vn ∈ Φi and vm ∈ Φj). Form a new particle from the mother particles, vnew = vn + vm. Eliminate particles n and m from their ensembles. If vnew = vi + vj ∈ [vk, vk+1], the new particle is allocated to Φk with probability Ps; otherwise, the new particle is allocated to subensemble Φk+1. The product probability parameter, Ps, is calculated using the mass-conserving equation for this landing interval: vkPs + vk+1(1 − Ps) = vnew. For example, in Figure 4, the next event was selected to be an aggregation of the two particles marked in black. These particles will disappear from the system to form a new particle that will be allocated either to subensemble 4 with probability Ps or to subensemble 5. The propensity or number of aggregation events per unit time is given by

(2)

⎛ 1 ⎞ R agg(i , j) = ⎜1 − δij⎟NNq i j (vi , vj)/ Vsim ⎝ 2 ⎠

where α is a user-specified constraint (typically 0.1−0.3).

3. FAST MC SIMULATION INSIDE EACH COMPARTMENT A second building block in the MCM algorithm is an efficient and accurate algorithm to model the particulate processes inside each compartment. Inside each compartment, the system is considered homogeneous and can be modeled by any existing MC method used for homogeneous processes. The method developed in refs 3 and 4 can decrease the computational cost by orders of magnitude without sacrificing accuracy. The PEMC and τ-PEMC simulations are event-driven and time-driven MC methods, respectively. The ideas in the construction of this method are discussed in refs 3 and 4. These algorithms are based

(3)

A nucleation event, Enuc, consists of inserting a new particle with volume vo in the first subensemble Φ1 (the first subensemble is set to represent particles of volume vo). The propensity function for this nucleation event is given by R nuc(En) = VsimJ(vo)

(4)

where J is the nuclei addition. 3.2. PEMC: Event-Driven Simulation of the PERP Model. The exact event-driven simulation of the PERP model consists of the following steps: Let Ei be the list of all possible events (Eagg(i,j), i = 1, M, j ≥ i, and Enuc). Let T be the total number of possible events (Ei, i = 1, ..., T). In this method, one event is 15487



Article

executed at a time. The time and the next event to be executed are selected stochastically using T

τ = −ln(r1)/∑ R(Ei)

(5)

i=1

The next event, Ef, is selected using the inverse method, where f satisfies f

T

f +1

∑ R(Ei) < r2 ∑ R(Ei) ≤ ∑ R(Ei) i=1

i=1

i=1

(6)

Set time to t = t + τ. 3.3. τ-PEMC: Time-Driven Simulation of the PERP Model. The approximated time-driven simulation of the PERP model consists of the following steps. In this method all events are fired at the same time many times. The time step is selected a priori. Given the same time list of possible events Ei, i = 1, ..., T, at each preselected time interval, τ, the number of times each event is fired is a Poisson random variable, Ks, with distribution PPD(Ks;REi,τ), where e − aτ PPD(K ; a , τ ) = (aτ )K K!

Figure 5. Arbitrary network of interconnected compartments. Each compartment is discretized in subensembles.

the number of particles to leave subensemble k is PBDout (Kout k,i ,pi,Nk,c(i)), where Kk,i is the number of particles in subensemble k that will leave compartment i. These particles are relocated to the temporary accumulation compartment and final compartment in the same k subensemble in each compartment. The modification of algorithm 1 is trivial (it basically consists of an outer adding loop over the subensembles). 4.1. Time-Driven Compartmental Simulation. The timedriven algorithm 3 is a time-split strategy that allows interfacing the PBF and the τ-PEMC method in a seamless way. A global time step control is implemented as described in this section. This will eliminate the complication of empty compartments. This procedure is described below. Algorithm 3: time-driven compartmental simulation 1. Initialize the particles in each compartment and set the time to zero. 2. Determine the time step Δt. 3. Apply τ-PEMC (section 3) to each compartment i using the time interval Δt. 4. Apply PBF(Δt) for each k subensemble using Nk,c(i) as the total number of particles for compartment i in step 2 of algorithm 1 and making all particle movements to the same k subensemble of the new compartments (accumulation or next compartment). 5. Increment the time t = t + Δt. 6. If t < tf, return to step 2; otherwise, end the simulation. Notice that the special capability of PBF to correctly model the particle transfer independent of the particle concentration is exploited in step 4. To facilitate the code implementation, the PBF is applied to each subensemble separately in the compartment (instead of the entire number of particles in the compartment). All of these subensembles have a large difference in the number of particles. The PBF method will handle these dilute subensembles with the same accuracy as the concentrated subensembles, thus eliminating any possible bias due to this partition in subensembles. 4.1.1. Determination of the Time Step. Let us define the compartment total propensity:

(7)

After all events have been fired (all possible aggregation events and nucleation events), the time is advanced to t = t + τ. The simulation is continued until a final time is reached. The accuracy and speed of the method depend on the selection of time τ during the simulation. A method for selecting τ was proposed in ref 25. Another variation of this method is to replace the Poisson distribution with a binomial distribution.26,27 In the τ-PEMC method, a binomial distribution is utilized and the time step is selected using a simple relation: T

τ = f /∑ R(Ei) i=1

(8)

where f is a parameter. Better time steps can be achieved using the formulas described in ref 25 at the expense of more complexity and computational cost.

4. MC SIMULATION OF COMPARTMENTAL POPULATION BALANCE MODELS Two CMC algorithms are developed to model particulate processes under heterogeneous conditions. The first algorithm is a time-driven algorithm, which is the main objective of this work (section 4.1). The second algorithm is an event-driven algorithm (section 4.2). This algorithm is discussed for the special situation of very dilute systems that are inherently stochastic.5 To implement the CMC algorithms, each compartment is discretized in subensembles (using the same discretization) to apply the PEMC or τ-PEMC algorithm. This is shown schematically in Figure 5 for the case of a three-compartment model. The subensemble k in compartment i is indexed Φk,c(i) with Nk,c(i) particles. The accumulation compartment used in the PBF algorithm also has the same discretization (also shown in Figure 3). To simplify the coding of the movement of particles between compartments, the PBF method is applied to each subensemble separately instead of to the total number of particles in the compartment. Applying the PBF to a subset of particles in a compartment is possible because the probabilities are given in terms of the flow statistics and time step, and not in terms of the particle concentration. In step 2 of the PBF algorithm,

R c(i) =

∑ R(Es ,c(i)) Es ,c(i)

15488

(9)



Article

where Es,c(i) is the list of all s events in compartment i and the total propensity is R tot =

∑ R c(i)

1. 2. 3. 4.

(10)

i

5.

The time step is selected using the formula

ΔtMC = f /R tot

(11)

This is eq 8 using the propensity of the total systems. For accurate simulations this time step should be bounded by the upper bounds imposed on the fluid flow: Δt PBF = min(UB h , UBf )

(12)

The actual time step is the minimum of these two quantities: Δt = min(ΔtMC,ΔtPBF). 4.1.2. Time Splitting. Similar to finite difference schemes, to increase the accuracy of this time-split strategy, the PBF step can be applied in two Δt/2 half-steps, PBF(Δt/2), one before and one after application of the τ-PEMC method (over the full step Δt). 4.2. Event-Driven Compartmental Simulation (Multiscale Problems). When considering a very dilute system or system with rare events (very dilute subpopulation), the output of some properties of interest may be a fluctuating variable with a PDF. In these cases, the τ-leap conditions may not apply and the use of the τ-leap algorithm may no longer be valid. The only accurate solution for this case is the exact event-driven simulation.3,5 If the rate of particle transfer between compartments is orders of magnitude larger than the rate of MC events, then the exact event-driven simulation will be very slow. The simulation will spend most of the time sampling particle transfers from compartment to compartment. One way to accelerate the simulation speed is to exploit the properties of the PBF algorithm in a hybrid strategy similar to the one developed in ref 5. This is possible again due to the capability of the PBF to correctly capture the probability distribution function of a small fluctuating population. Without this property, the hybrid strategy cannot be implemented because the transfer process will generate errors in the number of particles belonging to the small population. These errors then propagate in the sampling of the time for the next event and the selection of the next event. To simulate the hybrid time-driven PBF/event-driven MC method, I consider a partition similar to the multiscale problem in ref 5. The “slow” events are the population balance model events inside each compartment. The dynamics of these events is modeled using the exact event-driven simulation. The “fast” events are the particle flows between compartments modeled using the PBF method. Since the fluid flow (“fast events”) changes the total propensity of the system, the time for the next event is therefore given by

∫t

6. 7. 8.

5. NUMERICAL TESTING OF THE PBF METHODOLOGY Although the validity of the method is proved by construction in Appendix A, the PBF method is tested by calculating the holdup time and number of visits to each compartment for three benchmark problems consisting of an arbitrary network of compartments. The PBF method to solve these problems is presented in Appendix B (algorithm B.1). The results are compared with the exact event-driven stochastic simulation developed in Appendix B (algorithm B.2). 5.1. Prototype Problem 1. Consider the three-compartment model in Figure 6. The flow rates are { f12 = 0.01, f13 = 0.05, f 23 =

Figure 6. Compartment network used in benchmark problems 1 and 2.

0.1, f 31 = 0.2}. The initial conditions consist of 20 000 particles in the first compartment {N1 = 20 000, N2 = 0, N3 = 0}. The simulation consisted on the average of 100 runs, and each run was terminated when the final time was reached, tf = 300. This problem was solved using algorithm B.1. These results are considered as benchmark (called the exact MC results). The same problem was solved using algorithm B.2 using different fixed time intervals. These results (labeled PBF(Δt)) are compared with the exact MC results. The average total number of visits, Mi, the average total holdup time, Ti, for exact MC, and the PBF(Δt) for Δt = 0.1, Δt = 1, and Δt = 3 are shown Tables 1 and 2. In Table 3, a second moment of Ti is provided. The deviation (%) reported in these tables is calculated using deviation(x) = 100|xPBF(Δt) − xexact MC|/ xexact MC. These results demonstrate an excellent agreement with very low deviation even when very coarse times are used. The probability of violating the hopping constraint was calculated using algorithm 2. As shown in Figure 7a, for the coarser time the hopping constraint is violated 5% of the time on average.

t+τ

R tot dt + ln(r1) = 0 0

(13)

To find the event to be fired at τ, first compartment f where the event occurred is found by solving f

f +1

∑ Rc(i) < r2R tot ≤ ∑ R c(i) i=1

i=1

Set the system initial conditions. Initiate or reset the residual equation: RES = ln(r1). Select the time step, Δt, for the fast particle flow. Integrate the residual jump equation: RES = RES + Rtot(t)Δt. If RES ≥ 0: a. Adjust the Δt parameter to match the time for the next slow event: Δt = −RES/Rtot(t). b. Find compartment j where next event occurs using eq 14. c. Find the actual event inside compartment j by solving eq 6 using the list of events Es,c(j), s = 1, ..., T. d. t = t + Δt (if t ≥ tf, then end the simulation). e. Return to step 2. Perform PBF(Δt) applied to each subensemble. t = t + Δt (if t ≥ tf, then end the simulation). Go to step 4.

(14)

where r1 and r2 are random numbers. After compartment f is found, eq 6 is solved to find the actual particle event to be fired. Algorithm 4: exact event-driven algorithm (slow or dilute population dynamics and very fast f lows) 15489



Article

Table 1. Ti for the Exact Simulation and the PBF Method at Different Times method

deviation (%)

compartment

exact

Δt = 0.1

Δt = 1.0

Δt = 3.0

Δt = 0.1

Δt = 1.0

Δt = 3.0

1 2 3

215.85 20.85 63.30

215.79 20.87 63.34

215.83 20.85 63.32

215.83 20.86 63.31

0.03 0.12 0.07

0.01 0.02 0.03

0.01 0.06 0.02

Table 2. Mi for the Exact Simulation and the PBF Method at Different Times method

deviation (%)

compartment

exact

Δt = 0.1

Δt = 1.0

Δt = 3.0

Δt = 0.1

Δt = 1.0

Δt = 3.0

1 2 3

13.66 2.16 12.88

13.66 2.16 12.88

13.66 2.16 12.88

13.66 2.16 12.87

0.02 0.05 0.02

0.01 0.04 0.01

0.02 0.02 0.01

Table 3. Second Moment of Ti method

deviation (%)

compartment

exact

Δt = 0.1

Δt = 1.0

Δt = 3.0

Δt = 0.1

Δt = 1.0

Δt = 3.0

1 2 3

47222.90 773.10 4378.26

47188.50 772.63 4378.92

47154.30 754.80 4331.91

47038.30 719.22 4233.94

0.07 0.06 0.02

0.15 2.37 1.06

0.39 6.97 3.30

Table 4. Ti for the Exact Simulation and the PBF Method at Different Times compartment

exact

Δt = 1.0

dev (%)

1 2 3

74.77 37.13 2888.10

74.71 37.23 2888.06

0.08 0.27 0.00

Table 5. Mi for the Exact Simulation and the PBF Method at Different Times compartment

exact

Δt = 1.0

dev (%)

1 2 3

29.90 7.46 29.87

29.90 7.45 29.87

0.00 0.12 0.01

Table 6. Second Moment of Ti Figure 7. Percentage of particles hopping more than one compartment (violating the hopping constraint): (a) network used in benchmark problem 1, (b) network used in benchmark problem 2, (c) network used in benchmark problem 3.

5.2. Prototype Problem 2. Consider the three-compartment model in Figure 6, but this time the flow rates are { f12 = 0.1, f13 = 0.3, f 23 = 0.2, f 31 = 0.01}. The initial conditions consist of 20 particles in the first compartment {N1 = 20, N2 = 0, N3 = 0}. The simulation consisted on the average of 2000 runs, and each run was terminated when the userspecified final time was reached, t = 3000. This is a very challenging problem because the very small number of particles in the system will induce an intrinsic statistical noise. In this case, the average residence time is a probabilistic variable with a probability distribution function determined by running the statistics of many simulations (2000 simulations). Furthermore, any change in particles will be considered a very large change where traditional leaping conditions are violated. This shows that this algorithm works regardless of a low number of particles and drastic changes in compartment numbers of particles, as long as the hopping constraint is satisfied. The results shown in Tables 4−6 show that, even for this extreme case, very accurate results are obtained. The distributions of the

compartment

exact

Δt = 1.0

dev (%)

1 2 3

5.94 × 103 1.74 × 103 8.34 × 106

5.85 × 103 1.71 × 103 8.34 × 106

1.46 1.67 0.00

average holdup times for each compartment are shown in Figures 8−10. These distributions are virtually indistinguishable. These results were analyzed using analysis of variance (ANOVA), confirming that both simulation results are statistically identical. 5.3. Prototype Problem 3. In this problem, a network with a larger number of compartments is studied. The problem consists of 40 compartments. For the first 39 compartments, the flow rates are f i,i+1 = 0.1, and for the last compartment, the flow rate is f40,1 = 0.1. The definition of the network is completed with the following flow rates: f10,3 = 0.05, f15,35 = 0.05, f 38,25 = 0.05. The initial conditions consist of 20 000 particles in the first compartment {N1 = 20 000, N2 = N3, ... = N40 = 0}. The final time is tf = 1000. To generate the statistics, 2000 simulations were performed. The average holdup time is shown in Figure 11. The step changes in holdup times are a result of the flow bypasses ( f10,3, f15,35, f 38,25). The PBF was able to capture these step changes very accurately. The deviation in average holdup time and average 15490



Article

Figure 8. Histogram of the average holdup time in compartment 1: (a) exact method, (b) PBF method.



needed for the implementation of the random execution procedure described in Appendix A).

number of visits for each compartment is less than 0.02%. The two solution methods were compared using ANOVA, confirming that both solutions (exact MC and PBF) are statistically equivalent. The three prototype problems were also solved using a larger number of carriers (Li = 2Ni), where half of the carriers were empty. The simulation results were identical to the case where all carriers are real particles (Li = Ni). Notice that in practice Li ≈ Ni even when enough empty carriers are added to accommodate the creation of new particles in MC steps (this extra allocation is only

6. NANOPARTICLE FORMATION UNDER HETEROGENEOUS CONDITIONS USING A TWO-COMPARTMENT MODEL In this section, the formation of metal nanoparticles under heterogeneous conditions is modeled using a two-compartment model introduced in ref 28. In this model, based on CFD simulations, the system is divided into two compartments. 15491



Article

Figure 11. Average hold time in a 40-compartment network problem.

One compartment with volume Vimp represents the small region around the impeller characterized by large energy dissipation rates, εimp. The second compartment represents the larger circulation region with volume Vcir and a small dissipation rate, εcir. These parameters are determined from detailed CFD simulations. These quantities are related to the bulk properties as follows: εV = εimpVimp + εcirVcir

(15)

V = Vimp + Vcir

(16)

in particular the formation of oversized particles. To model the formation of secondary particles, a simplified model is implemented in which there is an initial concentration, cp,0, of primary particles in the reactor (very fast addition and nucleation). The primary particles grow by aggregation using the Brownian and turbulent kernels: q(vi , vj) =

The exchange between compartments is characterized by a volumetric flow, Q. The two-compartment model is shown in Figure 12. These two zones (a very small zone with a high energy

3kBT 1/3 (vi + vj1/3)(vi−1/3 + vj−1/3) 3vf ε 1/3 (vi + vj1/3)3 + 1.29 νf

(14)

A more complete model that includes a more complete kernel and autocatalytic production of primary particles is discussed in ref 33. This model can easily be incorporated in this scheme as well. In the simulations the primary particle radius is set to 10 nm. The initial concentration of primary particles in the reactor is cp,0 = 1.0 × 1019 M−3. The simulation volume is set: Vsim = 104/cp,0. The kinematic viscosity is ν = 1.0 × 10−6 M2/s. The values of the simulation parameters are summarized in Table 7. Table 7. Compartment Parameters: Vimp = χ1Vsim, Vcir = χ2Vsim, εimp = χ2ε, εcir = χ4ε RPM

χ1

χ2

χ3

χ4

ε

f 21

20 100 400

0.37 0.08 0.04

0.63 0.92 0.96

2.25 7.35 12.6

0.26 0.44 0.54

1.25 × 10−4 0.016 1.0

0.09 0.14 0.20

The particle size distribution (PSD) for the three cases is shown in Figure 13 (average of five simulations). As the number

Figure 12. Two-compartment model.

dissipation rate and a large zone with low energy dissipation) have been identified experimentally and computationally.29,30 Another extension to this model is a three-compartment model to include the effect of baffles as well.31 The two-compartment model is utilized to study the effect of disparate values in dissipation energy inside the reactor during the formation of metal nanoparticles. The main formation mechanism of metallic spherical particles is the production of primary particles that then aggregate to form secondary particles.32,33 On the basis of this mechanism, the formation of spherical metallic nanoparticles (silver or gold) is modeled using an aggregation population balance equation. Depending on the chemistry, the final particle size distribution and morphology are very sensitive to the process type and conditions.34 The goal is to understand the role of the high-intensity zone in the aggregation of primary particles to form secondary particles,

Figure 13. Effect of turbulence on the final PSD.

of impeller revolutions per minute (RPM) is increased, the PSD becomes broader. To understand the importance of the 15492



Article

This type of analysis can be extended to include more detailed zonal models than account for the effect of micro- and mesomixings and the injections of feeds. This analysis can also be applied to study other systems, such as the formation of quantum dots.35 These extensions are matters of future research.

high-energy zone in this broadening, the following parameters are recorded during the MC simulation: The percentage of total collisions that occur in compartment 1, Call 1 (%), the percentage of total collisions that produce particles larger than a cutoff radius, max Cr>r 1+2 (%), where rmax is the cutoff radius, and the percentage of total collisions that produce particles larger than a cutoff radius r>r and occur only in compartment 1, C1 max (%). These types of parameters are possible to calculate because of the capability of MC methods to follow individual events and particle history. These values are summarized in Table 8 for the three cases.

7. CONCLUSIONS Efficient and robust compartmental Monte Carlo methods for time-driven and event-driven simulations were developed to study population balance models under heterogeneous conditions. This type of analysis is a very important tool during scaleup and optimization of industrial-scale systems. In this work the PBF was introduced, which is a critical building block for the CMC methods. This is a time-driven method based on the construction of a binomial process of particles inside a CSTR flow configuration. Theoretical analysis and numerical simulations prove that this is a very robust and efficient method to transfer particles between compartments. Different from other leaping strategies, this method is valid for all concentrations, including compartments with a very dilute population of particles and independent of the complexity of the network. Furthermore, the PBF method captures the stochastic details of the exact stochastic simulation of particles moving between compartments. The capability of the PEMC and τ-PEMC methods to perform very fast simulations makes the overall CMC algorithms very efficient. The time-driven CMC method was utilized to study the effect of large variations in dissipation energy inside the reactor on the particle size distribution. A two-compartment model with parameters determined using CFD simulations was utilized. The aggregation of primary particles model was utilized with and without the influence of colloidal stability. It was found that turbulence induces broadening of the particle size distribution but most of the dynamics occur away from the impeller. It may be concluded that, since the area of very high intensity near the impeller has a secondary impact on primary particle aggregation, the main driver is the bulk energy input to the system. It was also found that a good stability factor can suppress the broadening induced by turbulence. These models will be extended in future work to include the effect of meso- and micromixing on the formation of primary particles.

Table 8. Statistics on the Collision History of the High-Dissipation-Energy Zone Q max Cr>r 1+2

max Cr>r 1

(%)

(%)

max Cr>r 1+2 (%)

max Cr>r (%) 1

15.78 5.75 5.33

0.04 0.17 34.20

0.01 0.03 5.01

14.40 5.35 4.39

0.00 0.00 2.57 × 10−3

0.00 0.00 8.57 × 10−4

(%) case 1 case 2 case 3

Q/W

Call 1

(%)

Call 1

As seen in this table, low values of Call 1 indicate that most of the max dynamics occurs away from the impeller. Comparing Cr>r with 1 r>rmax C1+2 also indicates that most of the large particles are formed away from the impeller as well. Therefore, the broadening is dominated by the increase of the average energy dissipation, ε, and not by the large variation between εimp and εcir. The same cases are studied but now under the influence of a colloidal stability using DLVO theory:

qW (vi , vj) =

qW (vi , vj) W (ri , rj)

(15)

The parameters used in refs 33 and 34 are used to calculate the stability factor, W(ri,rj). The only parameter modified was the Debye−Huckel parameter, κ, where a value of 3 nm was utilized. The effect of turbulence on PSD broadening is reduced as shown in Figure 14. The statistics in Table 8 show that the collision

■

APPENDIX A. CONSTRUCTION OF THE STOCHASTIC PBF MODEL I start the analysis considering the simplest scenario possible, a well-mixed compartment of volume V in isolation with N nonreacting particles, and with a steady volumetric flow in and out, q. Let f be the dimensionless volumetric flow rate q/V. Furthermore, the inlet stream of this idealized system carries no particles. Now consider a given time interval of interest, Δt. The number of particles leaving the compartment during Δt is a random variable with a binomial distribution. This is because a binomial process consists of a possible number of indistinguishable outcomes out of a maximum number of possible outcomes (given the probability for a single outcome to happen). In this case, the outcome is that the particle leaves the compartment, the maximum number of outcomes is the number of initial particles in the compartment, and the probability of an outcome happening is the probability of a particle leaving the compartment in a period Δt, which is p = fΔt. Therefore, the number of particles leaving the compartment during Δt is a random variable, Kout, which is a sample from a binomial distribution (eq 1) with

Figure 14. Effect of the stability factor on the final PSD.

history of the high-dissipation-energy compartment is very similar to that of the simulations without stabilization (Call 1 , %), but the number of collisions that produce large particles is minimized in both zones. As the number of rpm keeps increasing, that effect of turbulence in making large particles is increased as well. From these simulations, it can be concluded that if good stability (dispersing agents, etc.) is added to the particles, broadening of PSD due to turbulence can be suppressed during scale-up. It can also be concluded that most of the events that lead to large particles occur away from the impeller. 15493



Article

Figure 15. Binomial process for a network that satisfies the hopping constraint.

interval Δt. These changes will change the maximum number of particles, Nc(i), that can be transferred during this interval, affecting the binomial process just constructed since the assumption of an “isolated” binomial process is no longer valid. This problem is resolved by considering a fixed number of “carrier” particles, Li ≥ Nc(i), where all carrier particles n ∈ {1, ..., Nc(i)} correspond to an actual particle originally in the system and the extra carrier particles n ∈ {Nc(i) + 1, ..., Li} are “empty” carriers that contain no real particles. If a new particle is formed during Δt, then the new particle is assigned to one of the empty carrier particles and the carrier becomes a filled carrier representing the real particle. If a real particle is destroyed, then the corresponding carrier becomes an empty carrier particle with no real particles assigned to it. If Li is large enough to accommodate these changes, then the particle dynamics induced by the ongoing MC steps is hidden in the partition between filled and empty carrier particles. In this case, the carrier particles leaving compartment i in a time interval Δt again become an exact binomial process with maximum number of outcomes Li and parameter pi: PBD(Kout i ;pi,Li) (independent of the underlying MC process). The concept of carrier particles is indicated in Figure 16. In most cases, the effect of the MC process can be

parameters p = fΔt and maximum number of outcomes N: PBD(Kout;p,N). Now consider a network of interconnected well-mixed compartments numbered 1, 2, ..., Ncomp with volumes Vi. The network description is completed with the specification of the steady flow rates between compartments, qij. Let f ij = qij/Vi and f i = ∑j f ij. Initially, in each compartment i, there are Nc(i) particles. In this case, the number of particles leaving each compartment in a given time Δt is not guaranteed to be a binomial process as in the isolated compartment case. To recover a binomial process for this generalized case, I impose a constraint on the analysis of the time interval Δt of interest such that the probability of a particle leaping more than one compartment in this time interval is zero or below a tolerance level, P≥2 < εh, where P≥2 is the probability of a particle visiting more than one compartment and εh is a userspecified tolerance. I call this constraint the particle hopping constraint. If I consider a time interval that satisfies the particle hopping constraint, then the number of particles leaving any compartment i is again a binomial process as explained next. The hopping constraint guarantees that new upstream particles entering compartment i during Δt do not leave i during this interval. Therefore, compartment i will behave in the same way as the isolated compartment where the new entering particles do not interfere with the process and can be ignored. Thus, the number of particles that can leave the compartment is still Ni(t), and the probability of an outcome (leaving i) is still pi = f iΔt (the two conditions that define a binomial process). This is best illustrated using the drawing in Figure 15 where there are three compartments in series. Focusing on compartment 2, there are Nc(2) = 7 potential particles (marked as circles) that can leave this compartment in the time interval Δt. However, during the same period there are three particles entering from compartment 1 into compartment 2, marked as crosses. Since Δt satisfies the hopping constraint, the probability of each cross moving to compartment 3 is zero, and thus, I can still consider the binomial process of the particles originally in compartment 2 (circles) with seven particles (ignoring the effect of the upstream particles (crosses)). Thus, p2 = f 2Δt and Nc(2)(t) = 7. Now consider a new complication where there is an intracompartment Monte Carlo process such that the number of particles in each compartment may be changing during the time

Figure 16. Concept of carrier particles. A filled carrier particle represents a real particle.

ignored since the time steps are usually selected such that they only produce a very small disturbance to the system. Therefore, 15494



Article

setting Li = Nc(i) is a very accurate approximation for most practical purposes and is the one implemented in this work. If there are different categories of particles, each category can be considered as a separate binomial process. Also, for the generic case, the flows can be time-dependent. It is assumed that Δt is small enough that the fluid flow can be considered constant in this interval. Therefore, this method can be applied to any type of fluid flow, steady or unsteady. This completes the construction of the PBF model.

■

6. If a final time is reached, end the simulation; otherwise, go to step 2.

APPENDIX C. τ-LEAP METHODS FOR HETEROGENEOUS SYSTEMS Although there are existing methods that use Poisson statistics (τ-leap methods) to transfer molecule species between compartments, they do not possess the properties of the PBF. The main reason is that these methods are based on a different set of assumptions. Furthermore, they have been developed to model a different type of problem. In particular, spatially heterogeneous τ-leap methods have been developed to model the diffusion reaction equations (∂c/∂t = D∇2c + R) to simulate a biological system with spatial concentration gradients.22,23 These methods, in their different variants, are based on adapting the standard τ-leap method (developed for homogeneous chemical reactions) to include the diffusion process. This is accomplished by approximating the diffusion process with a set of first-order chemical reactions of the form Ci → Cj, where Ci is the concentration of a given chemical species in compartment i. The conversion from Ci to Cj in these reactions represents the transfer of the chemical species from compartment i to compartment j. The rate of this artificial reaction is given in terms of the diffusion coefficient and the grid dimensions (k = D/l2). These equations are added to the list of chemical reactions, and then the standard τ-leap algorithm is applied to the extended set of reactions. The τ-leap method is based on the assumption that all reaction rates (including diffusion rates) remain constant in the selected time interval. This is achieved when the time step is small enough such that the concentration in the compartment remains nearly constant after the firings of the reactions during the time interval. This is the so-called τ-leap condition.24 This condition is satisfied when the concentration in the compartment is high enough such that the change in concentration due to diffusion and reaction is a very small perturbation when compared with the concentration before the time step. Therefore, due to this constraint, the accurate implementation of this method requires a strategy for time step selection and a set of control parameters to switch to exact event-driven simulation when the number of molecules is below a threshold.22,23 These methods cannot be applied to very dilute systems. The PBF eliminates this fundamental difficulty. This is because the PBF is constructed to model fluid flow (instead of diffusion) and it is based on imposing a different type of constraint called the hopping constraint (instead of the τ-leap condition). The ability of the time-driven PBF to accurately simulate the transfer of particles regardless of the concentration of particles is a key property that differentiates this algorithm from the τ-leap algorithms. This property has important implications in modeling of small populations or subpopulations. It is also instrumental in the interface between the τ-PEMC method and the PBF.

■

APPENDIX B. EXACT AND PBF METHODS TO GENERATE HOLDUP TIMES An exact event-driven stochastic simulation is presented, exploiting the simple fact that the time for exiting the wellmixed compartment is a stochastic variable with an exponential distribution, f(τ) = exp( f iτ). The algorithm follows the history over a fixed time period of all particles, N, in the compartmental network. The problem is specified with the compartment network flows, the number of particles, and the initial number of particles in each compartment. The algorithm is utilized to compute the average holdup time in each compartment i, Ti, and the average number of visits to compartment i, Mi. The procedure is summarized in the following algorithm. Algorithm B.1: Event-driven stochastic simulation of residence time distribution 1. For each particle, n = 1, ..., N: a. Initialize Ti = 0, Mi = 0 for all compartments i. b. Place a particle in its initial compartment i. Mi = 1. c. Determine the exit time, τi, from τi = −ln(r)/f i, where r is a uniform random number in the interval [0, 1]. d. Update the time t = t + τi. e. Record the accumulated holdup time for the current compartment: Ti = Ti + τi. f. Pick another random number, r. Find the next compartment j, such that j is the smallest integer j satisfying rf i ≥ ∑ k=1 f ik. g. Record the visit to the new compartment j: Mj = Mj + 1. h. If a final time is reached, go to the next test particle; otherwise, go to step 2. 2. End loop over all particles. 3. Calculate the average holdup time and the average number of visits for each compartment i: Ti = Ti /N, Mi = Mi/N, where N is the total number of particles in the system. Traditionally, the discrete time Markov chain method is used to generate the RTD (residence time distribution) for a given network of continuous systems.20,21 In this method, a small constant time step is selected and the simulation is carried out at constant time interval increments. The method proposed in algorithm 5 does not depend on the discretization time since it is an exact simulation procedure that samples from the exponential distribution of each compartment. To test the PBF algorithm in terms of its ability to generate correct holdup times and numbers of visits to each compartment, the following procedure is implemented. Algorithm B.2: holdup times using PBF(Δt) 1. Initialize the system: set Ti = 0 and Mi = 0. 2. For each particle in compartment i, set Ti = Ti + Δt. 3. Execute PBF(Δt). 4. Redistribute: for every particle entering compartment j count the visit, Mj = Mj + 1. 5. Update the time: t = t + Δt.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel: (919) 248-5246. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Garcia, A. L.; van den Broek, C.; Aertsens, A.; Serneels, R. A Monte Carlo method of coagulation. Physica 1987, 143A, 535−546.

15495



Article

(2) Liffman, K. A direct simulation Monte-Carlo method for cluster coagulation. J. Comput. Phys. 1992, 100, 116−127. (3) Irizarry, R. Fast Monte Carlo methodology for multivariate particulate systemsI: Point ensemble Monte Carlo. Chem. Eng. Sci. 2007, 63, 7649−7664. (4) Irizarry, R. Fast Monte Carlo methodology for multivariate particulate systemsII: τ-PEMC. Chem. Eng. Sci. 2007, 63, 7665−7675. (5) Irizarry, R. Stochastic simulation of population balance models with disparate time scales: Hybrid strategies. Chem. Eng. Sci. 2011, 66, 4059−4069. (6) Zhao, H.; Kruis, F. E.; Zheng, C. Reducing statistical noise and extending the size spectrum by applying weighted simulation particles in Monte Carlo simulation of coagulation. Aerosol Sci. Technol. 2009, 43, 781−793. (7) Kruis, F. E.; Wei, J; Zwaag, T.; Haep, S. Computational fluid dynamics based stochastic aerosol modeling: Combination of a cellbased weighted random walk method and a constant-number MonteCarlo method for aerosol dynamics. Chem. Eng. Sci. 2012, 70, 109−120. (8) Cheng, J.; Yang, C.; Mao, Z. S. CFD-PBE simulation of premixed continuous precipitation incorporating nucleation, growth and aggregation in a stirred tank with multi-class method. Chem. Eng. Sci. 2012, 68, 469−480. (9) Alexopoulos, A. H.; Maggioris, D.; Kiparissides, C. CFD analysis of turbulence non-homogeneity in mixing vessels, a two compartment model. Chem. Eng. Sci. 2002, 57, 1735−1752. (10) Wells, G. J.; Ray, W. H. Methodology for modeling detailed imperfect mixing effects in complex reactors. AIChE J. 2005, 51, 1508− 1520. (11) Guha, D.; Dudukovic, M. P.; Ramachandran, P. A.; Mehta, S.; Alvare, J. CFD-based compartmental modeling of single phase stirred. AIChE J. 2006, 52, 1836−1846. (12) Moullec, Y. L.; Gentric, C.; Potier, O.; Leclerc, J. P. Comparison of systemic, compartmental and CFD modelling approaches: Application to the simulation of a biological reactor of wastewater treatment. Chem. Eng. Sci. 2010, 65, 343−350. (13) Irizarry-Rivera, R.; Seider, W. D. Model-predictive control of the Czochralski crystallization process. Part I. Conduction-dominated melt. J. Cryst. Growth 1997, 178, 593−611. (14) Irizarry-Rivera, R.; Seider, W. D. Model-predictive control of the Czochralski crystallization process. Part II. Reduced-order convection model. J. Cryst. Growth 1997, 178, 612−633. (15) Zheng, Y.; Johnston, D.; Berwick, J.; Chen, D.; Billings, S.; Mayhew, J. A three-compartment model of the hemodynamics response and oxygen delivery to brain. NeuroImage 2005, 28, 925−939. (16) Ninawe, P. R.; Hatziavramidis, D.; Parulekar, S. J. Delivery of drug macromolecules from thermally responsive gel implants to the posterior eye. Chem. Eng. Sci. 2010, 65, 5170−5177. (17) Zauner, R.; Jones, A. G. On the influence of mixing on crystal precipitation processesApplication of the segregated feed model. Chem. Eng. Sci. 2002, 57, 821−831. (18) Díez, M. D.; Fjeld, M.; Andersen, E.; Lie, B. Validation of a compartmental population balance model of an industrial leaching process: The Silgrain® process. Chem. Eng. Sci. 2006, 61, 229−245. (19) Portillo, P. M.; Muzzio, F. J.; Ierapetritou, M. G. Hydrid DEMcompartment modeling approach for granular mixing. AIChE J. 2007, 53, 119−128. (20) Pippel, W.; Philipp, G. Utilization of Markov chains for simulation of dynamic chemical systems. Chem. Eng. Sci. 1977, 32, 543−549. (21) Berthiaux, H.; Mizonov, V. Applications of Markov chains in particulate process engineering: A review. Can. J. Chem. Eng. 2004, 82, 1143−1168. (22) Marquez-Lago, T. T.; Burrage, K. Binomial tau-leap spatial stochastic simulation algorithm for applications in chemical kinetics. J. Chem. Phys. 2007, 127, 104101:1−9. (23) Iyengar, K. A.; Harris, L. A.; Clancy, P. Accurate implementation of leaping in space: The spatial partitioned-leaping algorithm. J. Chem. Phys. 2010, 132, 094101−1−14. (24) Gillespie, D. T. The chemical Langevin equation. J. Chem. Phys. 2000, 113, 297−306.

(25) Gillespie, D. T.; Petzold, L. R. Improved leap-size selection for accelerated stochastic simulation. J. Chem. Phys. 2003, 119, 8229−8234. (26) Tian, T.; Burrage, K. Binomial leap methods for simulating stochastic chemical kinetics. J. Chem. Phys. 2004, 121, 10356−10356. (27) Chatterjee, A.; Vlachos, D. G.; Katsoulakis, M. A. Binomial distribution based τ-leap accelerated stochastic simulation. J. Chem. Phys. 2005, 122, 024112. (28) Alexopoulos, A. H.; Maggioris, D.; Kiparissides, C. CFD analysis of turbulence non-homogeneity in mixing vessels: A two-compartment model. Chem. Eng. Sci. 2002, 57, 1735−1752. (29) Baldyga, J.; Bourne, J. R. Calculation of micromixing in inhomogeneous stirred tank reactors. Chem. Eng. Res. Des. 1988, 66, 33−38. (30) Bakker, A.; Myers, K. J.; Ward, R. W.; Lee, C. K. The laminar and turbulent flow pattern of pitched blade turbine. Trans. IChemE 1996, 74A, 485−491. (31) Vakili, M. H.; Esfahany, M. N. CFD analysis of turbulence in a baffled stirred tank, a three-compartment model. Chem. Eng. Sci. 2009, 64, 351−362. (32) Park, J.; Privman, V.; Matijevic, E. Model of formation of monodispersed colloids. J. Phys. Chem. B 2001, 105, 11630−11635. (33) Irizarry, R. Simulated dynamic optical response strategy for model identification of metal colloid synthesis. Ind. Eng. Chem. Res. 2010, 49, 5588−5602. (34) Irizarry, R; Burwell, L.; Leon-Velazquez, M. S. Preparation and formation mechanism of silver particles with spherical open structures. Ind. Eng. Chem. Res. 2011, 50, 8023−8033. (35) Leon-Velazquez, M. S.; Irizarry, R.; Castro-Rosario, M. E. Nucleation and growth of silver sulfide nanoparticles. J. Phys. Chem. C 2010, 114, 5839−5849.

15496


Fast Compartmental Monte Carlo Simulation of Population Balance

Recommend Documents