Mathematical Development and Comparison of a Hybrid PBM-DEMDescription of a Continuous PowderMixing Process

is paper describes the development of a multidimensional population balancemodel (PBM) which can account for the dynamics of a continuous powder mixing/blending process. e PBM can incorporate the important design and process conditions and determine their effects on the various critical quality attributes (CQAs) accordingly. e important parameters considered in this study are blender dimensions and presence of noise in the inlet streams. e 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. e 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.


Introduction and Background
Although the pharmaceutical industries must satisfy strict production speci�cation norms imposed by regulatory authorities, mainly due to inefficient control strategies [1,2] and the nonpredictive effects of input parameters, the �nal products obtained are oen 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 �uid processing units due to the absence of set of governing equations derived from the �rst principles which can describe granular �ow under speci�c conditions. e interactions of the particles with surrounding particles, �uid, 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 oen 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 Re�neries, Petrochemicals and Food) have adapted state of the art simulation techniques and satisfy their �nal product requirements due to their continuous modes of production and well understood �ow 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 �ow 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 de�ned. Various soware 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 �nite number of particles are considered to interact via several contact and non-contact forces. e 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][23][24][25]. It was coupled with computational �uid dynamics for describing particle-�uid interactions [26,27] and continuum models [27]. Various rotational mixers [28][29][30], helical mixers [31,32] and rotor type mixers [33][34][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 particleparticle or particle-�uid interactions and dynamic behavior of granular materials. is paper addresses model development of a continuous mixing process using gPROMS (tm) as the platform. is 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.). e critical quality attributes have not been obtained directly from EDEM. e 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. ese 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] soware package which allows both steady state and dynamic simulation runs. e 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. e 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.

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][43][44][45] but not for mixing till date.

Population Balance Equation.
e 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. e term ( ) ( )( ) accounts for the rate of change of particle distribution due to change in an internal coordinate (e.g., particle size). e term ( ) ( )( ) accounts for the rate of change of particle distribution with respect to spatial coordinates. ℜ formation and ℜ depletion stand for particles forming and depleting, respectively. In this work, free-�owing 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 �ow 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 Journal of Powder Technology 3 the multi-scale coupling of PBM with DEM of free-�owing particles and the population balance equation is reduced to:

Multidimensional Population Balance Equation for
Continuous Mixing. A multi-dimensional formulation of the PBE is considered. erefore, 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. e 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. e terms and represent the velocities in axial and radial directions, respectively. In�ow and out�ow terms are added accordingly to account for particles entering and exiting the mixer. In�ow is the rate at which the components are fed to the system. It is a constant value over time. If there are max × max compartments then out�ow can be represented as ∑ = ∑ max = ( max ) , where is the forward axial velocity.

Modeling Technique and Model
Outputs. e 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 �usti�ed in literature as well �17]. Particles are treated as discrete entities and their exchange between any two compartments is simulated. e 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. e 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 �xed axial location at any given point of time. Hence the mass balance of a single component can be simpli�ed according to e 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.
is 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.
e mixing performance is de�ned 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). ese CQAs should be regulated and controlled in order to achieve the desired mixing efficiency. ese 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. e 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 max is 2. max and max stand for the maximum number of grids in axial and radial direction, respectively. e homogeneity of samples retrieved from the out�ow is measured by calculating the variability in the concentration. e relative standard deviation (RSD) of tracer concentration measures the degree of homogeneity of the mixture and is given as represents total number of compartments ( = max × max ). i is the index to represent the compartment. is the concentration of component A at any compartment i. avg 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 �ow. 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 �nd the RTD: (i) the �ow in the blender is well mixed� (ii) the powder elements entering the blender simultaneously �ow with constant velocity and leave the blender at same time.

Numerical
Technique. e formulated PBM is a multidimensional hyperbolic partial differential equation (PDE). e PDE was �rst discretized using a central �nite 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.

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

DEM Simulation.
In EDEM (tm), a commercial blender (Gericke GCM250 (tm)) with impeller blades in alternating forward and backward orientation was simulated. e 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). e 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. e simulation was run for 260 seconds. e 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. ese data were then used to obtain the RSD as a function of blender length and time, rate of out�ow, component A composition at discharge, and RTD. In order to calculate the net out�ow, a bin was formed at the discharge. e 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 Journal of Powder Technology 5 time step. A code was written to �nd the out�ow in terms of total number of particles being discharged per time step. e 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. us the total rate of particle �ow 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 �owrates 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 10 × 10 bins both axially and radially. e 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. e 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 �atten 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. Aer the steady state was achieved, the particles that were fed to the blender within a 1-second window were tagged. ese particles were tracked until they crossed the weir at the outlet of the mixer. e time taken by each tagged particle to cross the weir was recorded as the residence time of the particle. e 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. e RTD, , was calculated by normalizing the area under the concentration-time curve.

gPROMs-Based PBM Simulation.
In gPROMs, the domain was discretised into 10 bins each in axial and radial coordinate axes. e width of the bin along the axial and radial coordinates were 33 mm and 10 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. e DEM output values (axial and radial particle velocity) have been extracted aer every 5 seconds (starting from 0 till the �nal time point 0 s) 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 1 0 s. Similarly the axial velocities were obtained. e 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, �rst component B stream was allowed to run through the blender. A pulse input of component A particles was introduced at 0 seconds aer the steady state was reached.

Hybrid PBM-DEM Model.
Combining DEM with PBM requires detailed understanding of both the models and the establishment of a well de�ned interface between them. e 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).

Effect of the Blender Dimensions.
Knowledge of minimum blender length and diameter required to ful�ll 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. e RSD decreases with time as well as over the axial length of the blender. e 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. is 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 �nal 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 �xed. 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].

Effect of Noise.
is section investigates how a mixer will respond to a possible perturbation in input �owrates. e usual source of disturbance at the inlet of the mixture is re�lling of the feeder [48]. e �owrate �uctuations 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 �owrate of the inlet stream which selects a value over a normal distribution. e 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.  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.
is 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.

Summary.
e 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. Figure 7 depicts the RTDs obtained from the DEM and gPROMs model. e 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. e 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. e 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. e RSD at discharge decreases over time for both the models as the system turns two segregated streams of components into more uniform blend. is 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. e steady state value is 0.5 since same amount of both the components were taken at the inlet. e fractional composition values for DEM again seem to highly �uctuate about the steady state whereas they change gradually in case of the hybrid PBM-DEM model. Figure 10 represents how the particle �owrate at the discharge vary with time. Both the gPROMs and DEM results �uctuate. is is because powder �ow cannot be explained on the basis of �rst principles unlike �uid �ow. In a continuum phase such as a �uid, a perfect steady-state is possible because the inlet and outlet �owrates can match exactly, which is not possible with discrete particles. Hence the concept of perfect steady state is not realised in powder system. ere 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 e�ect on the powder �ow. But the �uctuations 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.

Comparison of Hybrid PBM-DEM with Full DEM Simulation.
e 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. is clearly demonstrated the efficacy of the PBM model as a tool for design, control and optimization of continuous mixing processes.

Comparison of Simulation Time between DEM and PBM.
e full DEM took 6.5 days running on a 4 core and 2 threads/core processor with a total of 8 workers. e 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 signi�cantly more, taking up to 90% of available RAM compared to the PBM which uses up 50% RAM. is 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. en the PBM can be modi�ed 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].

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. e 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. e ultimate aim is to effectively use the PBM for control and optimization of blending.