A Multi-Scale Hybrid CFD-DEM-PBM Description of a Fluid-Bed Granulation Process

: In this study, a hybrid multi-scale model has been developed for a continuous ﬂuid bed wet granulation process by dynamically coupling computational ﬂuid dynamics (CFD) with a discrete element model (DEM) and population balance model (PBM). In this process, the granules are formed by spraying the liquid binder on the ﬂuidized powder bed. The ﬂuid ﬂow ﬁeld has been solved implementing CFD principles and the behavior of the solid particles has been modeled using DEM techniques whereas the change in particle size has been quantiﬁed with the help of PBM. The liquid binder droplets have been modeled implicitly in DEM. A detailed understanding of the process aids in the development of better design, optimization and control strategies. The model predicts the evolution of important process variables ( i.e. , average particle diameter, particle size distribution (PSD) and particle liquid content) over time, which have qualitative similarity with experimentally observed trends. The advantage of incorporating the multi-scale approach is that the model can be used to study the distributions of collision frequencies, particle velocity and particle liquid content in different sections of the ﬂuid bed granulator (FBG), in a more mechanistic manner


Introduction, Motivation and Objectives
Granulation is widely applicable in industries which deal with powder handling processes (e.g., food, pharmaceutical, catalyst, fertilizer industries etc.).During granulation process, fine powder particles form agglomerates which are granules of larger size with improved properties (e.g., flowability, uniform composition etc.).In pharmaceutical industries, wet granulation is an important unit operation in the downstream tablet manufacturing process, since it has a significant effect on the mechanical properties of the tabletting material (e.g., hardness, dissolution rate).As a result, the understanding of the process dynamics of this particular unit operation is critical.An inefficient operation of the granulation process leads to high batch rejection rates (if operated in batch mode) or high recycle ratio (if operated in continuous mode) [1].Since it is necessary to closely monitor the product quality at every step in powder handling processes, the Food and Drug Administration (FDA) has introduced the principles of Quality by Design (QbD) and Process analytical technology (PAT).These guidelines allow consistent building of product quality at every step of the manufacturing framework and help reduce production of inferior quality products.As the pharmaceutical industries are transitioning towards implementation of QbD and PAT guidelines, in order to improve the processing efficiency, there is a need to develop more science-based models which can be used to capture detailed process dynamics.Such a process model can be used for virtual experimentation and optimization, prior to testing in the actual plant.A mechanistic model developed from the first-principles can be more efficient in capturing the dynamics of a process compared to empirical or statistical models [2].Since fluid bed granulation process consists of both solid and fluid phases, it is desired to model both the phases individually and understand how they interact with each other.The fluid phase can be treated as a continuum and the flow dynamics can be described with the help of continuity equation and equation of motion.On the other hand, solid particles are discrete entities which will require the implementation of DEM techniques in order to capture their flow dynamics.Such types of flow which include both fluid and solid phases can be simulated efficiently with the help of CFD-DEM coupling.CFD will account for fluid flow whereas the same can be done for the solid phase using DEM.Establishing a connection between the flow field of the fluidizing medium and the contact pattern of the particles which in turn will dictate the aggregation rate will provide a detailed process dynamics of the fluid bed granulation process.A multi-scale modeling scheme in this case is required.The advantage of developing a multi-scale model is that it stores the information from various scales (i.e., micro, meso and macro) and couples the information to simulate the process dynamics.Continuum models are associated with macro scale simulation, volume averaged equation over a cell with grid size larger than the particle size is an example of meso scale simulation and in micro scale simulations, the flow of each and every individual particle is considered [3].For the first time, this work introduces a framework which uses macro scale information from a CFD model (based on continuity equations), meso scale information based on a PBM and particle level information from DEM to simulate the process dynamics.

Objectives
The purpose of this study is to develop a hybrid CFD-DEM-PBM model for fluid bed granulation process with the following objectives: • Present a hybrid CFD-DEM-PBM framework using dynamic two-way coupling.
• Incorporate multi-scale information such that the model can be used to study the detailed process dynamics.• Study the heterogenous particle velocity distribution and liquid binder distribution.
• Study the evolution of average particle diameter and particle liquid content with time.
The CFD model has been solved using Fluent TM (version 14.5.0)(ANSYS Inc., Canonsburg, PA, USA), DEM has been solved using EDEM TM (version 2.5.1)(DEM Solutions, DEM Solutions Ltd.Edinburgh, UK) and in order to link the PBM with DEM, a user-defined library model has been built using the application programming interface within EDEM TM .

Background
In order to facilitate agglomeration during granulation process, a liquid binder is added and granulation takes place as a result of three rate processes (1.wetting and nucleation; 2. consolidation and aggregation and 3. breakage and attrition) [4].Based on these rate processes, a granule can be classified by three internal properties (e.g., size, liquid content and porosity).In a FBG, the particles are fluidized in order to ensure maximum contact between the liquid and the powder particles.At first the particles are rendered wet when they come in contact with the liquid binder and form granule nuclei.As the wet particles collide with one another, they form liquid bridges resulting in aggregation of small sized particles to form a larger sized particle.Similarly, particle breakage occurs because of the presence of shear stresses, compressive and tensile forces within the granulator due to particle-particle interaction or particle-wall (vessel interior wall) interaction, which can break a particle into smaller fragments.
The particle wetting rates in a FBG is highly dependent on the particle flow pattern which in turn is affected by the flow field of the fluidizing gas [5].The flow behavior of the powder within the granulator also depends on the geometry of the vessel.Several attempts are being made to understand the granular flow pattern but many aspects of it (i.e., a more detailed study of the vessel geometry, fluid flow field, collision frequencies etc.) using sophisticated modeling tools and techniques are yet to be developed [6].The physical properties of the granules that are important in deciding the compactibility of the granules and which also affect the mechanical strength of the tablets, are particle liquid content, particle size, particle surface area and pore diameter.Porosity relates to the void spaces present in the particle and it is an important particle property because the surface liquid content of a particle will depend on its porosity and the surface liquid content will further decide the particle's aggregation rate.Porous particles can easily deform during compaction to create new bonding surfaces, thus affecting the tablet hardness [7].Hardness of the tablet also depends on the granule surface area [8].An appropriate design and operation of the granulation process can help to control the PSD and also help to achieve the desired flow characteristics of the particles.Improved granule characteristics will lead to efficient operation of the other unit operations present in the tablet manufacturing framework [9] and also help to achieve better product quality.
PBMs have been used extensively in order to model granulation process [10,11].PBMs are used to calculate the rate processes (e.g., aggregation, consolidation and breakage) during granulation, but certain information, used in calculation of mechanistic aggregation and breakage kernels (i.e., effect of particle properties, spatial effect, collision frequency, particle velocity etc.) cannot be determined from a PBM alone.A DEM on the other hand is able to determine these entities.Such a PBM-DEM coupled model for particle aggregation has been reported [12].A one way coupling consisting of PBM-DEM hybrid model has been reported [13,14] in case of continuous mixing, in which the PBM consists of spatial discretization only (as there is no size change taking place during mixing).
Since a CFD-DEM approach has been used to study several complex particle-fluid flow systems [15][16][17], it is an effective tool to model fluid bed granulation.Interested readers are referred to the review presented by [18] on the discrete particle modeling of fluid beds.A combined CFD-DEM approach was first presented by [19] for simulating plug flow through horizontal pipes.Since then, there have been several works reported in scientific literature [20][21][22][23].Liquid-solid interactions have been studied in case of gas bubble formation in a gas-liquid-solid system [24][25][26].Similar studies have been reported in case of gas and liquid fluidized beds in the bubbling regimes where a rigorous two way coupling framework has been introduced to explain the fluid-particle interaction [27,28].Fernandez et al. used CFD-DEM coupling in case of centrifugal separation systems [29].Previous studies have been carried out in CFD-DEM coupling in case of granulation and are well documented [5,30].A direct numerical simulation model for a three phase flow in case of wet granulation, consisting of solid, liquid and gas has been developed by [3].The liquid-gas flow has been solved using one of the CFD multi-phase solvers, and the solid particle flow has been solved by DEM.
A CFD-PBM framework has been already implemented in order to model size change taking place in a particle-fluid system [31,32] and a PBM derived from kinetic theory of granular flow (KTGF) has been developed by [33] for fluidized melt granulation.Hydrodynamic modeling of wet granulation process in a FBG using CFD principles has been reported by [34].However there is a major scope still open to study a CFD-DEM-PBM coupled framework, which will capture a detailed interaction between the phases and quantify the size change.Such a multi-scale model formulation has been reported by [35] for agglomerate breakage in fluid bed.The present work aims at capturing the particle aggregation during a granulation process in a FBG.The framework can be further developed in future by adding particle breakage along with the aggregation.

Multi-Scale Model Development
This section illustrates the mathematical equations and modeling techniques of each domain of the coupled framework.The CFD model calculates the fluid flow-field and the DEM adds it to the force-field calculation for each particle.The PBM obtains the various information required to calculate the rate processes from DEM to quantify the size change and update the PSD.

CFD Model for the Fluidizing Medium
The fluid phase has been simulated using FLUENT TM .The flow model chosen is that of a laminar viscous flow.The governing equations for incompressible flow are given in Equations ( 1) and (2) [5,36]: Navier-Stokes Equation and Equation of continuity where is the solid volume fraction, u is the fluid velocity, ρ is the fluid density, p is flow pressure, µ is the viscosity, S p is the source term and g is gravity.The fluidizing medium in this case has been considered to be air and the fluid velocity at the inlet has been kept constant at 30 m/s.The fluid velocity at the inlet has been scaled as a function of the particle size in DEM to simulate observed experimental fluidization in a granulator.The inlet velocity has been determined by running several simulation trials such that the particles are fluidized and the effect of gravity can be nullified.In real granulation process, the particle size is in microns.But in the simulation, the initial particle size is in millimeters (i.e., 2 mm diameter), which is much higher than the initial particle size used in real granulation process.As the simulation bed weighs more, therefore a higher air velocity is required to counter the effect of gravity.The solution methods used for calculating the momentum and particle volume fraction is first order upwind [37].A first order Implicit scheme has been selected for the transient formulation.The boundary conditions have been set as follows: • Flow near wall is laminar and the velocity varies linearly with the distance from wall.
• A no slip boundary condition has been set at the wall.
• A velocity inlet boundary condition has been used for the air entering the geometry.
• An outlet-vent boundary condition has been used at the geometry exit.

Discrete Element Model
DEM essentially uses Newton's laws of motion (as shown in Equations ( 3) and ( 4)) to simulate the particle force fields (namely contact forces and body forces [38]).The contact forces are due to particle-particle or particle-boundary (vessel internal wall) contacts whereas the body forces are any external force fields (i.e., gravity or fluid flow field in case of fluid bed granulation) acting on the particles.The net force F total is calculated for each particle at a time interval which is approximately in the order of 10 −5 or 10 −6 s and the new particle state is calculated by numerically solving Newton's law and Euler's equations of rotational motion.The simulation uses a damped Hertzian normal contact model with Mindlin-Deresiewicz/Coulomb friction tangential force model.A detailed discussion on the contact models along with the governing equations have been provided by [38].
The DEM has been simulated using EDEM TM .The geometry of the FBG has been imported within EDEM TM and the initial PSD has been created.The initial PSD is uniform such that all particles are of fixed size (1 mm in radius).A virtual particle factory plate has been introduced within the geometry in the lower half and the particles are created at a generation rate of 50,000 particles per second.The liquid addition has been captured in EDEM TM .The liquid droplets have been modeled as solid particles with similar properties as the initial powder particles.Liquid particles have been created continuously at a rate of 50,000 particles per second from a virtual plate placed in the upper half of the granulator.It should be noted that the liquid addition has been started at time equal to 0.2 s in order to let the particles fluidize first.The feed rate and the number of particles have been set in a way such that the simulation speed is fast and a reasonable granulation is achieved within the simulation time (when compared qualitatively with the experimental results [39,40]).In a FLUENT TM -EDEM TM coupled simulation, it is important that the particle size is comparable with the mesh size (i.e., the particle size shouldn't be very large compared to the mesh size) [41].For example in this work the primary particle diameter is 2 mm and the mesh size is 4 mm.These parameters have been fixed by running several trial simulations.
The change in particle size and new PSD is calculated using the PBM in the subsequent time steps as will be explained later (in Section 3.3).Whenever a liquid particle collides with a powder particle, it is deleted from the system and the liquid content of the particle increases by the liquid particle volume.If two liquid particles collide with each other, a single droplet with volume equal to the total volume of the two particles is formed.It should be noted that as the liquid particles are deleted from the system immediately upon contact, therefore EDEM TM does not consider the contact forces due to these particles in the simulation.The EDEM TM tracks the collisions occurring in between different particles and the information is stored in form of an array, which is made use by the PBM for calculating the new population distribution function at every population balance time step.Table 1 gives the material properties and the other parameters for particle-particle and particle-wall interaction.The density of the powder material is 1030 kg/m 3 (which is similar to that of a pharmaceutical active ingredient (S)-Ibuprofen [42,43]).Detailed information on the exact range of parameter values (particle-particle and particle-wall interaction parameters of EDEM TM ) for S-Ibuprofen is not available in literature, hence the values have been adapted from [38], who simulated fines in a continuous mixer, and from the work of [44].

Population Balance Model for FBG
The Population Balance Equation (PBE) can be expressed as shown in Equation ( 5) [45]: Here, F(x,z,t) is the population distribution function, x is the vector of internal co-ordinates used to express the particle size, z is the vector of external co-ordinates used to represent spatial position of the particles and t is the time.The term ∂ ∂x F (x, z, t) dx dt accounts for the rate of change of particle distribution due to change in particle size.The term ∂ ∂z F (x, z, t) dz dt accounts for the rate of change of particle distribution with respect to spatial co-ordinates.R f ormation and R depletion stand for particles being formed and depleted respectively due to aggregation and breakage.
The present study considers the particle size change only, therefore the internal coordinates (represented using solid and liquid volume) have been retained and the spatial coordinates have been dropped from the PBE.At present only the particle aggregation has been accounted for and breakage has been neglected.Equation (5) has been modified accordingly to obtain Equation (6) as shown below: The first and second terms of the equation represent evolution of the particle distribution function with respect to time and rate of liquid addition respectively.s and l are the volumes of solid and liquid (per particle), respectively and dl dt is the liquid addition rate.The liquid addition has been captured implicitly in EDEM TM (as explained in the previous Section 3.2).Therefore the liquid addition term has been omitted from Equation ( 6) and the PBE has been further modified as shown in Equation ( 7): The aggregation rate process is defined in Equation ( 8) to Equation ( 10) where β(s, l, s , l , t) is the aggregation kernel, as given in Equation (11).The aggregation kernel has been defined as a function of the collision frequency (C) and collision efficiency ψ.The collision frequency (which is a function of particle size [46]) is calculated based on the number of collisions occurring between the particle groups which can be obtained from EDEM TM based on the particle properties.
β(s, l, s , l , t) = C(s, l, s , l , t)ψ(s, l, s , l ) The collision frequency is an important parameter.As previously mentioned, it is already known to be a function of particle size [46], but a study has been conducted by the authors showing collision frequency as a function of the PSD as well.A summary of this study has been provided as an appendix in this manuscript for the benefit of the readers.The effect of the particle size is captured in the collision frequency which in turn controls the aggregation rate.Therefore, the aggregation rate kernel of this model depends on the PSD at any point of time.Collision frequency can be calculated as shown in Equation ( 12) [46]: In the above equation, N coll is the total number of collision between two particle types represented by solid and liquid bins (F (s, l, t) and F (s , l , t)) during the time interval ∆t.
A simple expression (as shown in Equation ( 13)) has been adapted for collision efficiency based on the works of Biggs et al.
Here LC stands for liquid content of the particles and LC min is the minimum liquid content required for the particle coalescence.The above expression essentially means that if the liquid content of any two particles is greater than or equal to the minimum value (specified as 0.2 in this case), then the particles may aggregate to form a new particle upon collision depending on the value of the collision efficiency (which has been kept at a constant value of 0.01 in this study).The liquid content and collision efficiency values have been chosen such that the process variables show qualitative similarity in trend when compared to experimental results [40].
The PBM has been implemented by creating a custom contact model and custom factory for particle aggregation and liquid addition within the EDEM TM simulation.The initial particles (both liquid and solid) which are created within EDEM TM has a uniform solid and liquid volume respectively.Since a particle consists of both solid and liquid part, therefore the internal coordinate has been discretized linearly based on the particle's solid and liquid volume (also referred to as "bins").

Information Exchange in the Coupling Framework
The information exchanged over the coupled network along with the model assumptions have been summarized in this section.The model assumptions have been listed below: • The PBM considers aggregation only, breakage and consolidation has not been incorporated since FBG processes are low shear processes with reduced consolidation and breakage (similar approach has been followed by [48]).• A simple aggregation kernel has been formulated based on collision frequency and collision efficiency (adapted from [46]).• The collision efficiency in the aggregation kernel is size independent, non-mechanistic and conditional based on the liquid content of the powder particles (adapted from [47]).• Liquid addition has been captured in EDEM TM by creating particles which get deleted from the system upon contact.
Figure 1 shows a schematic of the main CFD-DEM-PBM coupling framework and the DEM-PBM framework has been magnified to show the ongoing steps within it.The CFD-DEM framework has been developed by coupling EDEM TM and FLUENT TM through the commercially available coupling interface between them.CFD simulates the flow-field of the fluidizing air.When the solution converges, the fluid flow-field is passed to the CFD-DEM coupling interface which then calculates the drag force acting on each particle.The calculated drag force is then transferred to the DEM solver which updates the particle flow-field to obtain the new particle positions (or state).It should be noted that the PBM has been coded within EDEM TM using the Application programming interface.The PBM time step is greater than the DEM time step.Therefore as the DEM solver performs the iteration, PBM waits till the PBM time step is reached.The number of collisions between each pair of solid and liquid bins over the time interval is recorded within the DEM simulation.This information is transferred to the PBM when the PBM time step is reached.The PBM uses this information to calculate the aggregation kernel as given by Equations ( 8)- (12).Each time the PBM is solved, a new PSD is calculated.This new PSD is then implemented within the DEM simulation for the subsequent DEM time step until the next PBM step is encountered.The DEM solver iterates until the CFD time step ends.The CFD-DEM coupling interface then takes the information of the new particle position (or state) from the DEM solver, updates the solid volume fraction in each fluid cell and passes the information to the CFD solver.The CFD solver again iterates over the next time step until the solution converges and the same steps are repeated.The EDEM TM time step should be always less than the FLUENT TM time step and it is suggested that the EDEM TM time step is kept between (1/10)th to (1/100)th of the FLUENT TM time step for stability [41].In this case, the CFD time step is 2.0 × 10 −3 s and the DEM time step is 3.25 × 10 −5 s.It means that the DEM solver iterates approximately 61 times within one CFD time step.The PBM time step is 0.25 s (i.e., PBM is solved every 0.25 s).The framework has been simulated for 2 s.Therefore the PBM has been solved 8 times in this simulation time.The PBM time step is an important factor and has been decided after running a few simulation trials.It is important to choose the time-step in a way such that: 1.A reasonable number of collisions occur among the particles between any two subsequent time steps.2. The PBM is solved a reasonable number of times such that there is a more consistent distribution of the particle size (as described in Section 4.3) If the PBM time step is too low, then the model will not be able to capture enough number of collisions which results in difficulty of computational tractability of PBM.A resolution study of the PBM time step and its effect on PSD has been carried out (details have been provided in the results and discussion section (Section 4.3)).The coupled framework has been run for a time period of 2 s which took approximately 6 h of CPU time, running on an Intel Core i7-2600 CPU processor (3.4 GHz) (Intel Corporation, Santa Clara, CA, USA) with 16 GB of RAM.

Model Outputs
This section presents a brief description of the model outputs which are PSD, average fractional liquid content and average diameter.
The PSD has been obtained by post-processing the simulation data.The diameter of each particle present in the system at any particular time point can be obtained from EDEM TM .Several size classes are defined by ranges of the diameter and each bin from the PBM are grouped in these size classes.The total mass of particles in each of these size classes is determined individually and each of them is further normalized by dividing with the overall mass of particles to obtain a mass frequency as seen in Equation ( 14).The mass frequency can be plotted with respect to the size classes in order to obtain the PSD.
where µ m is the mass frequency of mth size class (which is a function of particle diameter), ρ is the density, ns is the number of solid bins, nl is the number of liquid bins, V m is the particle volume in m th size class, nm is the number of size classes and D ij is the particle diameter corresponding to the ith solid bin and jth liquid bin.
Similarly the average diameter (D avg ) can be calculated as shown in Equation ( 15): The average fractional liquid content (x avg ) can be calculated as shown in Equation ( 16):

Results and Discussion
This section includes the discussion on the model geometry and the results obtained from the coupled framework.A brief section on setting up the coupled model has been included which enumerates the step wise set-up of the simulation.

Simulation Procedure
This section lists the procedure followed while setting up the coupled simulation, in a nutshell.
1.The geometry has been made using ANSYS Design Modeler TM .2. The geometry has been meshed using ICEM-CFD TM .3. The mesh file has been imported within FLUENT TM .4. The mesh has been converted into Polyhedra domain.5.The gravity is defined in the correct direction and a transient simulation is selected.6.The flow model has been selected to be viscous laminar.7. The coupling server has been started.8.The FLUENT TM is coupled with EDEM TM for the desired fluid domain by selecting the Eulerian-Eulerian option.9.The coupling server will automatically import the geometry with the specified direction of gravity in EDEM TM and set the source terms in x-momentum, y-momentum and z-momentum calculation.
The value of the simulation parameters of the coupling interface has been set as follows: • Sample points: The number of points used by FLUENT TM to calculate the volume fraction of the fluid cell.This value has been set at 10, which means that a large particle can transfer its volume between 10 cells.This particular parameter decides the stability and speed of the simulation.A higher value of sample point may increase the stability but decrease the simulation speed.• Relaxation factor: The relaxation factors again help with stability and convergence of the solution.Reducing the value helps to increase stability and achieve convergence.Both momentum-MTM-under-relaxation factor and volume under-relaxation factor have been set at 0.7.
10.The inlet fluid velocity has been defined as 30 m/s.
11.The custom contact model and custom factory (for PBM calculation) have been imported within EDEM TM .12. The material properties, particle-particle and particle-wall interaction parameters as given in Table 1 have been set in EDEM TM .13.The initial PSD has been created in EDEM TM .14.The liquid particles have been created in EDEM TM (the liquid addition starts at 0.2 s). 15.Once the EDEM TM simulation is set up, initialize the solution in FLUENT TM .16. Run the calculation.

Model Geometry
The measurement details of the model geometry has been given in Figure 2. Figure 2 also presents the mesh, which has been done in ICEM-CFD TM .The geometry has been meshed using the tetra/mixed meshing tools with a grid size of 4mm.The mesh has been converted into polyhedra domain after importing it in FLUENT TM .As shown in Figure 2, the geometry has been divided into five fluid domains or cell zones (inlet domain, fluid inlet domain, main fluid domain, fluid outlet domain and outlet domain).The CFD-DEM coupling has been performed for the cell zone (main fluid domain).Every time the simulation leaves Fluent TM and enters EDEM TM , all the relevant data in the linked zones are passed to EDEM TM .Reducing the number of linked zones will reduce the amount of data transferred between the two softwares and thus increase the simulation speed.Therefore the simulation has been linked for the "main fluid domain" only.Although the coupling has been done for the "main fluid domain", the particles are free to move around throughout the geometry, but only the data from the "main fluid domain" are passed to EDEM TM .The results/plots reported in the manuscript are for the "main fluid domain" only.It should be noted that the spatial co-ordinates have not been considered in the PBM.The main idea is to present a CFD-DEM-PBM framework which is able to capture the granulation dynamics.Therefore the assumption is that the "main fluid domain" is representative of the whole geometry.

Multi-Scale Model Results
The multi-scale model has been simulated to obtain the average diameter and liquid content of the particles as a function of time.This model can be also used to study the distribution of particle velocity and particle liquid content, which will help to understand the flow pattern and liquid content of the solid particles.Since the aggregation rate depends on the liquid content of the particles, therefore it is important to study these aspects.It will also help in understanding the impacts of the equipment design and improve it for a better operational efficiency.
In order to analyze the particle flow pattern and liquid distribution, the "main fluid domain" has been divided equally into two sections (upper half and lower half).Each half has been divided into several compartments and the DEM simulation has been post processed to obtain the particle velocity in each of these compartments.As shown in Figure 2, the geometry has been divided into 30 grids in both x and z direction and 1 grid in y direction, such that there are 30 × 1 × 30 compartments in total.Figure 2 gives an illustration of the compartments in 2D as seen from the top view.Figure 3 shows the contour plots of the particle velocities (i.e., the resultant velocity of the three components v x , v y and v z ) averaged over different time intervals.Figure 3a,b present the particle velocities for both lower and upper halves, averaged over time 0 s-0.5 s and 0.5 s-1 s respectively.It can be seen that the velocity of the particles increases with time.However the lower section shows a more uniform distribution of velocity compared to the upper section.The plot for upper section in Figure 3a shows a non-uniform distribution with high velocity present sporadically in few locations and the velocity values in most of the compartments are closer to zero.This is because the particles are being injected into the system from a virtual plate placed in the lower half and it takes about 0.4 s for the bed to be fully fluidized.Similarly Figure 3c,d present the particle velocity distribution in both the halves averaged over time 1 s-1.5 s and 1.5 s-2 s respectively.These plots also show that the particle velocity increases with time as the bed gets fully fluidized and presents a more uniform velocity distribution in the upper half.It can be also noted that an absolute steady state in the particle velocity is not being realized.This is due to the mechanics of the granulation process (as the partcles interact with each other and the boundary in numerous ways, it is difficult to obtain an absolute steady state).
Figure 4a,b show the particle liquid content distribution averaged over 0 s-0.5 s (initial time of simulation) and 1.5 s-2 s (final time of simulation).Figure 4a shows a non-uniform distribution over a higher range for the upper part, with the values in most of the compartments lying close to zero (similar observation has been made in the velocity distribution as well).The distribution becomes comparatively uniform with time as seen in Figure 4b.The liquid content of the particles in the upper part is more than the liquid content of the particles present in the lower part.This is because the liquid droplets are being injected into the system from a virtual plate located near the upper half of the granulator, hence only a few of them are able to penetrate through the void space and reach the lower half.This shows that the binder liquid distribution is highly dependant on the particle flow pattern, which in turn also affects the aggregation rate (particle liquid content is an important factor in deciding the granulation rate processes).

Liquid addition starts
As previously mentioned, the amount of liquid present in the particles is an important attribute as it controls the rate process, which in turn will decide the increase in particle size.The snapshots (Figure 5) of the liquid content have been taken at 1 s and 2 s.The particles have been color coded based on their fractional liquid content.It can be seen that more particles are getting wet with time.This is because the particle velocity increases with time and a more uniform velocity distribution is obtained throughout the granulator (also seen from the contour plots) resulting in more contacts between the liquid and solid particles.So the average liquid fraction (i.e., the ratio of liquid volume to the total volume of a particle, where total volume is equal to the summation of liquid volume and solid volume of a particle) increases with time and levels off gradually as can be seen in Figure 6 (plot of particle liquid content versus time).Since the particles are free to move around (in all domains), therefore there are particles exiting the "main fluid domain".Initially there is a steep increase in the liquid content because the bed takes some time to be fully fluidized.As a result the number of particles in the "main fluid domain" is still increasing and more particles are getting wet with time.But as the bed is fully fluidized, an approximate steady state is being realized in the "main fluid domain".It should be noted that in this model, the liquid droplets which have been created are of comparable size to the solid particles.Therefore an assumption has been made that one liquid particle can wet only one solid particle.Under certain circumstances, where the powder particles are much smaller than the liquid particles, it is possible that one liquid droplet wets more than one powder particle (depending on the liquid to solid volume ratio).This feature will be included in the coupled framework in future.
Figure 7 presents the snapshot of the relative particle diameter taken at 1s and 2s.It can be seen that the number of large sized particles increase with time.This can be again explained on the basis of the particle movement pattern.With time, the velocity distribution becomes more uniform and the contact between the particles increases which results in aggregation.Figure 8 is a plot for the average diameter as a function of time.It is seen that the average diameter of the particles increase with time.Figure 9 is a plot for the PSD taken at different time points.It can be seen that a gradual progression towards larger sized particles is being realized during granulation.The trends observed in this study are qualitatively consistent with experimental results [39,40].These studies show that the granulation is being achieved as the particle liquid content increases.They present bimodal particle size distributions that broaden over time by starting with finer primary particles.
In order to show the effect of PBM time step, two more simulatons have been run with PBM time step equal to 0.2 s (PBM has been solved 10 times) and 0.4 s (PBM has been solved 5 times) and compared with the base case (PBM time step equal to 0.25 s, where PBM has been solved 8 times).Figure 10 presents the PSDs obtained at the final time point (time = 2 s) from each simulation.It can be seen that the PSD for the largest PBM time step (0.4 s) has a lower frequency of very large particles (due to the presence of the highest peak).This is because the PBM has been solved only 5 times during the simulation thus limiting the opportunities for large particles to collide and aggregate with each other.On the other hand, PSDs with smaller PBM time step (0.2 s and 0.25 s) show a more consistent distribution suggesting more accuracy.Therefore it is essential to perform a comparative analysis of the model predicted data obtained by running the simulation with different PBM time step, with experimental data.

Conclusions
This work presents a multi-scale hybrid model for granulation, using information from different scales of CFD, DEM and PBM techniques.Granulation has been extensively modeled following a PBM approach, which groups a lump of particles in different classes (based on size, porosity etc.), but this particular framework tracks the flow-field of each particle with the help of DEM.The CFD aspect helps to add in the drag force acting on the particles due to the flow field of the fluidizing medium.The aggregation kernel is a function of the collision efficiency and collision frequency, which in turn depends on the particle size distribution.An empirical expression, based on the liquid content has been used to determine the collision efficiency.Future effort will be to introduce a mechanistic approach in determining the collision efficiency based on relative velocity of the colliding particles, particle mass and liquid content.The PBM which has been implemented is two-dimensional which can track the change in particle size as well as liquid content of the particles.The liquid addition has been captured implicitly in the DEM and the aggregation kernel has been decided based on the information provided by the DEM on collision frequency and particle liquid content.PBM calculates the new PSD at every PBM time step, which has been implemented in the DEM by removing old particles and creating new ones.The model can be used to study the material flow pattern, equipment geometry (see if dead zones are present) and key performance criteria (i.e., average diameter and liquid content).The present framework considers particle aggregation only, but it will be further extended with the addition of breakage and consolidation processes in future.The model can be updated by addition of spatial co-ordinates as well.This model can be used to study the effect of material properties as well, which can be an interesting future investigation.This coupled model is detailed as it is able to store information from different scales and can be used as an effective tool to understand the process dynamics of a fluid bed granulation process.

Figure 1 .
Figure 1.Schematic of the coupled multi-scale framework.

Figure 2 .
Figure 2. Model Geometry and mesh of the FBG.

Figure 5 .Figure 6 .
Figure 5. EDEM TM snapshots of liquid content.(a) Particle liquid content at time = 1 s; (b) Particle liquid content at time = 2 s.

Figure 9 .
Figure 9. PSD at different time points.
] Mass frequency [−] PSD at PBM Time step=0.2 s PSD at PBM Time step=0.25 s PSD at PBM Time step=0.4 s

Figure A1 .
Figure A1.PSD used in DEM simulations to determine collision rates.