#### Abstract

This paper describes the development of a multidimensional population balance model (PBM) which can account for the dynamics of a continuous powder mixing/blending process. The PBM can incorporate the important design and process conditions and determine their effects on the various critical quality attributes (CQAs) accordingly. The important parameters considered in this study are blender dimensions and presence of noise in the inlet streams. The blender dynamics have been captured in terms of composition of the ingredients, (relative standard deviation) RSD, and (residence time distribution) RTD. PBM interacts with discrete element modeling (DEM) via one-way coupling which forms a basic framework for hybrid modeling. The results thus obtained have been compared against a full DEM simulation which is a more fundamental particle-level model that elucidates the dynamics of the mixing process. Results show good qualitative agreement which lends credence to the use of coupled PBM as an effective tool in control and optimization of mixing process due to its relatively fewer computational requirements compared to DEM.

#### 1. Introduction and Background

Although the pharmaceutical industries must satisfy strict production specification norms imposed by regulatory authorities, mainly due to inefficient control strategies [1, 2] and the nonpredictive effects of input parameters, the final products obtained are often nonuniform with a high level of variability with respect to product quality [3]. Moreover, the behavior of powder processing units are not well characterized as compared to the fluid processing units due to the absence of set of governing equations derived from the first principles which can describe granular flow under specific conditions. The interactions of the particles with surrounding particles, fluid, or equipment wall is quite complex to understand, model and manage. Bulk material behavior is decided by the interactions among individual particles at microscale, which is chaotic. Hence often the pharmaceutical industries have to follow a univariate trial and error approach for their process development. However efforts are being made in order to introduce science-based holistic development of process and product by using Quality by Design (QbD) and Process Analytical Technology (PAT) tools [4, 5].

Continuous manufacturing offers many advantages such as better process understanding and control. Several other chemical industries (e.g., Petroleum Refineries, Petrochemicals and Food) have adapted state of the art simulation techniques and satisfy their final product requirements due to their continuous modes of production and well understood flow processes [1]. A batch process with multipurpose equipment was shown to be more efficient than with single purpose equipment [6] because of more adaptibility of the equipment for various purposes. Continuous manufacturing processes are suitable for easy scale-up operations [2], improve product quality, and reduce the operating cost [7, 8]. Risks which are otherwise associated with solid handling and nonpredictive manufacturing also are mitigated [9]. It is important to understand each and every unit operation which will form a part of the continuous process in order to achieve the desired operational level. Common manufacturing steps which one can come across in a typical pharmaceutical industry are feeding, mixing, granulation, milling, tableting, and coating. Continuous processes can be designed and optimized with the help of flow sheet modeling, which has been shown to be a robust and detailed tool for simulating a real plant [10].

Mixing is one of the most important unit operations because the blend quality is primarily decided by this step. Mixing is brought about due to the particle velocities and velocity gradient within the blender when two or more distinct bulk material particles come into intimate contact. Segregation however can occur and induces variability in the mixture composition [11]. A model-based approach is a good way to understand the blender dynamics provided the parameters lie within the design range and are well defined. Various software packages and programming languages (e.g., ASPEN, gPROMS, DEM, MATLAB) are available to aid in this effort.

Among the several other modeling approaches which exist in literature, for example, Monte-Carlo methods [12], continuum and constitutive models [13], statistical models [14, 15], compartment models [16, 17], RTD models [18, 19] and hybrid models [20, 21], discrete element modeling (DEM) is one of the fundamental modeling approaches that is able to capture the particle level physics. In DEM, each particle is treated as a discrete entity where the trajectory of the particles is tracked and the collision between particles is modeled. A finite number of particles are considered to interact via several contact and non-contact forces. The translational and rotational motion of each particle follows Newton's Laws of Motion. DEM has been used as a tool for capturing the mixing dynamics by various kinds of systems [22–25]. It was coupled with computational fluid dynamics for describing particle-fluid interactions [26, 27] and continuum models [27]. Various rotational mixers [28–30], helical mixers [31, 32] and rotor type mixers [33–35] have been studied with the help of DEM. It has been used as an effective tool in order to study various particulate and multiphase systems but it still needs further developments in order to address the current research needs [27, 36] which would in turn help to understand the nature of particle-particle or particle-fluid interactions and dynamic behavior of granular materials.

This paper addresses model development of a continuous mixing process using gPROMS (tm) as the platform. This work builds upon a previous study by the authors [37] which followed a steady state approach to develop a population balance model (PBM) to describe the dynamics of the mixing process. In the present work a dynamic and novel hybrid mixing model which combines DEM with PBM has been compared with the data from a full DEM simulation run on EDEM (tm) (DEM Solutions Ltd.). The critical quality attributes have not been obtained directly from EDEM. The output of EDEM has been post processed in order to extract the CQAs. Moreover due to high computational time requirement, DEM is not an effective tool for control and optimization. Hence in this work an offline coupling between DEM and PBM has been considered. DEM simulation can give information regarding particle properties (particle velocity) which are on particle level. These particle level information can be fed to the PBM from which the macroscopic variables (RSD, RTD, blend composition, etc.) affecting the entire unit operation (mixing) can be extracted. In this way the model incorporates multiscale information and illustrates one-way coupling where DEM provides the velocity information and is combined with the PBM which simulates key blend attributes as a function of time and thus applies the microscopic properties from particle level in order to capture macroscopic properties which characterize the mixing performance of the entire blender.

gPROMs (tm) is a robust and fast equation-oriented [38, 39] software package which allows both steady state and dynamic simulation runs. The process model can be built by developing fundamental mathematical expressions relating various physical and chemical variables/parameters without specifying the order in which these equations need to be solved. The motivation behind this work is to present a more dynamic system which updates the particle properties at regular interval of time and generates the CQAs taking information from DEM, that can also be simulated in a realistic time period to facilitate design, control, and optimization.

#### 2. Mathematical Model Development

A population balance model (PBM) has been developed for capturing the dynamics of mixing. Population balance models have been used in case of other particulate handling processes such as crystallization [40, 41] and granulation [42–45] but not for mixing till date.

##### 2.1. Population Balance Equation

The generalised form of Population Balance Equation is given as [46] Here, is the population distribution function, is the vector of internal state variables on which the population distribution function depends on, and is the vector of external coordinates used to depict spatial position. The term accounts for the rate of change of particle distribution due to change in an internal coordinate (e.g., particle size). The term accounts for the rate of change of particle distribution with respect to spatial coordinates. and stand for particles forming and depleting, respectively. In this work, free-flowing particles have been considered and hence there will be no size change associated due to formation or depletion which would have been the case if the flow was cohesive thus potentially leading to particle aggregation. Formation and depletion terms may be added easily as per our previous study [37]. However, the focus of this study is the multi-scale coupling of PBM with DEM of free-flowing particles and the population balance equation is reduced to:

##### 2.2. Multidimensional Population Balance Equation for Continuous Mixing

A multi-dimensional formulation of the PBE is considered. Therefore, the equation for continuous mixing (see (3)) can be written as is the spatial coordinate in the axial direction and is the spatial coordinate in radial direction. The model deals with mixing of two components A and B such that represents component A and is for component B. is the particle number density which varies with spatial location inside the blender and type of particle. The terms and represent the velocities in axial and radial directions, respectively. Inflow and outflow terms are added accordingly to account for particles entering and exiting the mixer. Inflow is the rate at which the components are fed to the system. It is a constant value over time. If there are compartments then outflow can be represented as , where is the forward axial velocity.

##### 2.3. Modeling Technique and Model Outputs

The mixer has been divided into multiple zones both in the axial and radial directions. Mixing can occur in both axial and radial direction by convection and dispersion. In a continuous blender, mixing takes place when the particles are moved about by the motion of the blades with the dispersive component being negligibly small as compared to the convective one. Such assumptions have been justified in literature as well [17]. Particles are treated as discrete entities and their exchange between any two compartments is simulated. The space is discretized into several compartments and it is assumed that homogenous mixing occurs in each of the compartments. Particles move from compartment to compartment in both axial and radial directions and this is governed by the axial and radial velocities. The exchange of mass between the compartments has been represented as number of particles. Particles can either move forward to the compartment ahead of it or backward to the compartment behind it. On the other hand, radial mixing conserves the total number of particles at a fixed axial location at any given point of time. Hence the mass balance of a single component can be simplified according to The above equation can be written for component A and B. Here, refers to the forward velocity in the axial direction, refers to backward velocity in the axial direction, and refers to the radial velocity.

This model does not require any of the particle properties such as diameter, density, geometry, and so on as input provided the velocities can be measured either experimentally or from detailed numerical simulations such as DEM. Once velocity parameters are selected, the simulation can be used to provide information about the dynamics and the outcome of the process. It can be used to predict mixing performance in terms of a relative standard deviation (RSD) of sample along the axial length at end time point, RSD at discharge as a function of time, blend composition at discharge, and residence time distribution (RTD).

It can be noted that DEM will capture the particle properties in the velocity values. So if the particle size distribution is changed the velocity distribution is also going to change, which in turn will change the response of the PBM. Hence the various particle properties can be varied in order to determine how it will affect the CQAs.

The mixing performance is defined in terms of certain critical quality attributes such as relative standard deviation (RSD), composition of A, which can be active ingredient of interest for a given process (), and residence time distribution (RTD). These CQAs should be regulated and controlled in order to achieve the desired mixing efficiency. These parameters can be found as In the above equation, the numerator stands for the total number of particles of component A which come out of the last compartments at any point of time. The denominator represents total number of particles of both the components A and B coming out of the last compartments at any point of time. Since the model involves mixing of two components, so the value of is 2. and stand for the maximum number of grids in axial and radial direction, respectively.

The homogeneity of samples retrieved from the outflow is measured by calculating the variability in the concentration. The relative standard deviation (RSD) of tracer concentration measures the degree of homogeneity of the mixture and is given as
represents total number of compartments . *i* is the index to represent the compartment. is the concentration of component A at any compartment *i*. is a spatial average of component A concentration.

Residence time distribution is a measure of the time spent by the particles within the blender. In other words it captures the non-ideality associated with the flow. RTD can be found as
In the above equation, stands for concentration of component A in the outlet stream at any time *t*. It is important to make the following assumptions in order to find the RTD: (i) the flow in the blender is well mixed; (ii) the powder elements entering the blender simultaneously flow with constant velocity and leave the blender at same time.

##### 2.4. Numerical Technique

The formulated PBM is a multi-dimensional hyperbolic partial differential equation (PDE). The PDE was first discretized using a central finite difference scheme of order 6 which was followed by using an implicit backward differential formula (BDF) technique to integrate the system of ordinary differential equations (ODEs). Both the discretization and integration were performed using gPROMS (tm) in-built functions, which are state-of-the art schemes that ensures stability of the overall system and minimal numerical errors and numerical diffusion.

#### 3. Results and Discussion

All simulations were performed using a desktop computer with a 2.94 GHz Intel (Core i7) processor with 8 GB RAM.

##### 3.1. DEM Simulation

In EDEM (tm), a commercial blender (Gericke GCM250 (tm)) with impeller blades in alternating forward and backward orientation was simulated. The details regarding the mixer blade geometry has been elaborated by Dubey et al. [47]. Figure 1 gives a snapshot of the mixer geometry on EDEM (tm). The length and diameter of the mixer were 330 mm and 100 mm, respectively. Equal number of particles each of component A and B were introduced into the mixer using two feeders discharging particles on either side of the inlet. Table 1 gives the details about the particle properties, particle-particle, particle-blade, and particle-wall interaction parameters used in the simulation. A feed rate of 1990 particles per second and an impeller speed of 250 rpm were maintained. Normal particle size distribution with a mean radius of 1 mm with 5 standard deviation were used. The simulation was run for 260 seconds. The simulation was postprocessed to obtain the axial velocities, radial velocities, and particle IDs. In a DEM simulation, each particle is assigned a unique number known as the particle ID. These data were then used to obtain the RSD as a function of blender length and time, rate of outflow, component A composition at discharge, and RTD.

In order to calculate the net outflow, a bin was formed at the discharge. The IDs of the particles present in this bin were obtained at each time step for one time frame. Few particle IDs might get repeated between any two consecutive time step because some particles stay in the bin for more than one time step. A code was written to find the outflow in terms of total number of particles being discharged per time step. The code compared the particle IDs of every two consecutive time step, eliminated the particle IDs encountered in the previous time step, and increased the particle count by one whenever a new particle ID appeared. Thus the total rate of particle flow at the outlet was calculated. In the next step, the particle IDs were obtained separately for each particle of component A and B at the discharge so that the individual flowrates could be calculated. From this information, component A composition and RSD can be calculated for every time step with the help of (5) and (6), respectively. In order to determine the variation of RSD with the blender length, the blender was divided into bins both axially and radially. The individual numbers of particles for components A and B in each bin were obtained at the last time step. RSD values averaged over the blender length were calculated using (6).

In order to calculate RTD, the blender run was simulated until the mass hold-up in the blender reached a near-constant value, indicating that a steady state has been achieved. The hold-up starts at zero at the beginning of the simulation and will rise quasi-linearly in the beginning as the particles collect inside the blender. Once the particles start exiting the blender, the curve will flatten and it turns into a horizontal line when steady state is reached. At this point the number of particles entering the blender is about the same as those exiting. In the DEM simulation with the parameters shown in the paper, steady state was typically achieved at around 50s. After the steady state was achieved, the particles that were fed to the blender within a 1-second window were tagged. These particles were tracked until they crossed the weir at the outlet of the mixer. The time taken by each tagged particle to cross the weir was recorded as the residence time of the particle. The simulations were run at steady state for time intervals long enough so that at least 95 percent of the tagged particles were retrieved at the outlet. A histogram was created using 1-sec time bins and the number of particles in each time bin was plotted against time. The RTD, , was calculated by normalizing the area under the concentration-time curve.

##### 3.2. gPROMs-Based PBM Simulation

In gPROMs, the domain was discretised into 10 bins each in axial and radial coordinate axes. The width of the bin along the axial and radial coordinates were mm and mm, respectively. It is important to determine and suitably incorporate the velocity values (i.e., the axial and radial velocities) into the PBM. Each compartment has its own radial and axial velocity values. The DEM output values (axial and radial particle velocity) have been extracted after every 5 seconds (starting from till the final time point ) as excel sheets and then imported in the gPROMS model as foreign object. Figure 2 gives a snapshot of the radial velocity values as obtained from EDEM for different compartments at one of the time points . Similarly the axial velocities were obtained. The velocity values were updated every 5 seconds in the gPROMs model for a total time of 260 seconds. It should be noted that a more frequent update of velocities could be implemented (given that in DEM the velocities are calculated each time in the order of micro-seconds. However, this study focussed on the computational efficiency of the model without compromising on too much accuracy of results. Equal number of component A and B particles were introduced in the blender. Equations (5)-(6) were used to determine the critical quality attributes. In order to determine the RTD, first component B stream was allowed to run through the blender. A pulse input of component A particles was introduced at seconds after the steady state was reached.

###### 3.2.1. Hybrid PBM-DEM Model

Combining DEM with PBM requires detailed understanding of both the models and the establishment of a well defined interface between them. The model generated CQAs as explained earlier depend on the input parameter space. In the following sections, we investigate the robustness of the hybrid model by varying few of these parameters (i.e., dimensions of the blender and introduction of noise in the feed rate).

###### 3.2.2. Effect of the Blender Dimensions

Knowledge of minimum blender length and diameter required to fulfill the CQA requirements is essential from an equipment design point of view. Figures 3(a) and 3(b) represent the RSD versus time (at the mixer outlet) and RSD versus axial length (at the end time point), respectively for change in blender length while diameter is kept constant. The RSD decreases with time as well as over the axial length of the blender. The axial length has been represented in terms of compartment number (1 to 10). It can be seen from the graphs that RSD decreases with increase in the blender length. This is because as length increases, mixture is retained within the blender for a longer time thus giving it more time to get mixed. And the final product obtained is more uniform with reduced variability. Similarly Figures 4(a) and 4(b) show how RSD varies with change in diameter of the blender when the length is fixed. It is seen that the mixture variability increases with increase in the diameter. Similar results have been obtained in a previous study by the authors [37].

**(a)**

**(b)**

**(a)**

**(b)**

###### 3.2.3. Effect of Noise

This section investigates how a mixer will respond to a possible perturbation in input flowrates. The usual source of disturbance at the inlet of the mixture is refilling of the feeder [48]. The flowrate fluctuations at the blender inlet should be minimized so that the properties of the output stream from the blender are not affected. It has been shown that a continuous mixer can dampen out variability from the feeder [49]. Noise was added by adding a variance term to flowrate of the inlet stream which selects a value over a normal distribution. The standard deviation of the normal distribution has been varied in order to get an idea of the maximum allowable perturbation such that the output stream properties are not changed. Figure 5 shows how the fractional composition of component A varies at the outlet with time as the degree of perturbation changes. It can be seen that all the cases except the one where the standard deviation is 0.3 almost overlap. Figures 6(a) and 6(b) represent how the RSD changes with change in degree of perturbation. From Figure 6(a), it can be seen that the RSD deviates slightly for standard deviation of 0.3 whereas the rest overlap. Change in RSD with respect to time is not very evident as degree of perturbation changes because RSD approaches zero from a high initial value for all the cases.

**(a)**

**(b)**

This means that the developed model is robust and the mixer can eliminate any disturbance of small magnitude present in the inlet stream provided the degree of perturbation is within range.

###### 3.2.4. Summary

The hybrid PBM-DEM model demonstrated good qualitative agreement with experimental studies as well as full featured DEM simulations [15, 18, 50]. A further experimental validation of the PBM has been carried out by the authors [51] which is quantitative in nature.

##### 3.3. Comparison of Hybrid PBM-DEM with Full DEM Simulation

Figure 7 depicts the RTDs obtained from the DEM and gPROMs model. The plot shows that there is good qualitative agreement between the methods. Residence time can be increased by increasing the length and decreasing the blender speed [37]. Increasing the length or decreasing the blender speed will have considerable effect on other CQAs and cost. Hence it is crucial to optimize the blender performance as a function of processing conditions and formulation properties. An RTD study may be helpful in such type of process optimization.

Other CQAs were extracted from gPROMs and DEM as described in the previous two sections. The values were normalized and then compared against each other. Figures 8(a) and 8(b) show how the RSD varies with blender length and time, respectively. The overall RSD decreases over the blender length for both DEM simulations and the PBM-based gPROMs simulation. It can be noted that in DEM there are spikes occurring in RSD for compartments 2, 4, 7, and 8. On the other hand the decrease in RSD in case of gPROMs model is smooth. The RSD at discharge decreases over time for both the models as the system turns two segregated streams of components into more uniform blend. This shows that there is qualitative agreement between the two models as far as the blending dynamics is concerned. Both the plots show that the DEM results are very noisy and this is an inherent property of the simulation which assumes large sized particles due to which relatively small number of particles reside in the blender or exit the blender at any moment of time. On the other hand the hybrid PBM-DEM model shows very gradual variation of the properties. Figure 9 represents the component A concentration of the mixture at the outlet as a function of time. The steady state value is 0.5 since same amount of both the components were taken at the inlet. The fractional composition values for DEM again seem to highly fluctuate about the steady state whereas they change gradually in case of the hybrid PBM-DEM model.

**(a)**

**(b)**

Figure 10 represents how the particle flowrate at the discharge vary with time. Both the gPROMs and DEM results fluctuate. This is because powder flow cannot be explained on the basis of first principles unlike fluid flow. In a continuum phase such as a fluid, a perfect steady-state is possible because the inlet and outlet flowrates can match exactly, which is not possible with discrete particles. Hence the concept of perfect steady state is not realised in powder system. There are several ways in which any two particles can interact with each other as well as with the blender wall. Hence the particle-particle interactions and particle-wall interactions will have a pronounced effect on the powder flow. But the fluctuations seem to reduce over time. Overall, it can be seen that within acceptable error due to numerical noise, the PBM simulation in gPROMS is able to qualitatively capture the dynamics of the mixing process as demonstrated by a full DEM simulation. Moreover the hybrid PBM-DEM model is less noisy and more gradual while reporting the values of the CQAs.

The PBM-DEM model has been simulated for the same time interval as the DEM simulation and the computational requirements are reported in Section 3.4. This clearly demonstrated the efficacy of the PBM model as a tool for design, control and optimization of continuous mixing processes.

##### 3.4. Comparison of Simulation Time between DEM and PBM

The full DEM took 6.5 days running on a 4 core and 2 threads/core processor with a total of 8 workers. The PBM simulation on the other hand took 30 minutes running on a single core processor using 1 worker. Moreover, the memory occupancy of the DEM is significantly more, taking up to 90% of available RAM compared to the PBM which uses up 50% RAM. This clearly demonstrates the efficacy of using the PBM for control and optimization as opposed to the full DEM simulation which is not amenable to provide signal feedback given the time (of the order of days) it takes to perform a simulation. DEM simulation can be run only once in order to extract the particle level data (particle velocity), which can be fed to the PBM. Then the PBM can be modified with run as many times as required to extract the required macroscopic scale variables which affect the overall unit operation and thus make control and optimization easier because of its lower time requirement. It should be noted that the current PBM simulation takes 30 min in a serial simulation. Parallel simulation of PBMs using multicore CPU computing has shown to be efficient in further reducing the computational time of simulating a PBM thus enhancing its utility in control and optimization [52, 53].

#### 4. Conclusions

A hybrid framework of multi-dimensional population balance model (PBM) and discrete element modeling (DEM) was developed. PBM coupled with DEM forms a basis of one-way coupling. Variations in several design and process parameters such as number of compartments, blender dimension, and presence of disturbance in the inlet streams were considered in order to test the robustness of the hybrid model. PBM was shown to be an effective tool for tracking the blending dynamics in terms of the key properties such as blend composition and RSD. The results thus obtained from the hybrid framework were compared with the DEM. It gave a good quantitative agreement with the trends as seen in DEM. Future work will focus on two-way coupling of PBM and DEM and validation with experimental data. The ultimate aim is to effectively use the PBM for control and optimization of blending.

#### Nomenclature

: | Fractional composition of component A, dimensionless |

: | Average spatial composition of component A, moles/ |

: | Composition of component A in ith compartment, moles/ |

: | Composition of component A at any time t, moles/ |

: | Residence time distribution, dimensionless |

: | Particle density, particles/ |

: | Total number of compartments, |

: | Counter for number of components, |

: | Time, |

: | forward axial velocity, |

: | backward axial velocity, |

: | radial velocity, |

: | Spatial coordinate in axial direction, |

: | Spatial coordinate in radial direction, . |

#### Acknowledgment

This work is supported by the National Science Foundation Engineering Research Center on Structured Organic Particulate Systems, through Grant NSF-ECC 0540855.