CFD modelling of powder flow in a continuous horizontal mixer

This work presents a continuum model for simulating the flow of powders inside continuous horizontal mixers. The challenge is to adopt a reliable rheological model that allows simulating granular flows accurately. We selected the μ (I)-rheology model. First, we considered a set of granular collapse experiments, showing that the model can successfully reproduce these flows


Introduction
Flow and mixing of powders are the backbone of many industries, such as those related to pharmaceuticals, food, cosmetics, and cement [1].The consistent homogeneity of the mixtures plays a significant role in the finished product quality.The mixing process has recently shifted from traditional batch units towards continuous mixers [2].The pharmaceutical industry is a good example of the batch-to-continuous mixing shift, because this industry can be regarded as the reference for the standards of mixture homogeneity [3].The switch enables the continuous manufacturing of drug products (CMDP) [4].In CMDP, all operations (e.g., dispensing, blending, high-shear granulation, drying, tabletting) are integrated into a single production train without stopping and starting between units [5].CMDP can result in improved and more consistent product quality and a more flexible supply chain, reducing scale-up activities, manufacturing costs as well as environmental impact [5,6].Nevertheless, in CMDP processes controlling the homogeneity of the mixer outlet mixture is challenging [4].
To control the quality of the outlet mixture in continuous mixers, engineers must understand how powders flow and mix inside these systems.In particular, understanding powder flow is essential for designing, optimizing, operating and controlling continuous mixers.Knowing in detail the flow characteristics, such as velocity profiles, local shear rates, and mixer loading, facilitates the design of new, more effective mixers with appropriate geometries and impeller configurations.In addition, because the flow characteristics directly influence the quality and uniformity of the resulting mixture, investigating powder flow patterns in continuous mixers, including the identification and elimination of stagnant or dead zones, permits optimizing the fluid dynamics and preventing operational bottlenecks [7,8].
Many experimental works have investigated the effect of mixer design and of various parameters (e.g., powder properties and mixer operating conditions) on the performance of continuous horizontal mixers [9,10] using a variety of experimental techniques, such as visual tracking (VT), positron emission particle tracking (PEPT), particle image velocimetry (PIV), electrical capacitance tomography (ECT), computed tomography (CT) and magnetic resonance imaging (MRI) [11].However, to obtain local properties from within the mixer, such as velocity fields, these approaches can be expensive or difficult (if not highly impractical) to use [12].In contrast, mathematical modelling is more convenient for studying mixer performance, offering detailed insight into the complex flow patterns of the powders within the mixer.Costeffective and scalable, simulations provide a flexible alternative to experimental methods, overcoming their limitations, allowing for easy parameter modification, and yielding valuable information on pressure, velocity and concentration fields.Nevertheless, this insight is only as good as the models used to obtain it; thus, developers must ensure that the models capture enough physical details (but not much more than required by the goals of the simulations), that the values of the model parameters (in particular of the rheological material properties) are sufficiently accurate, and that the model predictions are reliable.Therefore, model validation is a necessary and key requirement.
At the most fundamental level, to predict the flow of powders, one should solve Newton's second law of motion for every particle in the powder [13].This approach, called the discrete element method (DEM), is accurate since the mean particle velocity and volume fraction fields can be calculated with minimum approximations [14][15][16][17][18].Moreover, since the granular medium is modelled as a collection of individual particles and not as a continuum, the model does not involve a stress tensor associated with the particle phase, and therefore no rheological problem of closure arises.This is a key advantage of the DEM modelling approach.
Blending processes in various types of mixers (e.g., rotating drums, bin blenders, twin-screw mixers, tote blenders, double-cone blenders, slant cone mixers, and ribbon mixers) have been successfully simulated using DEM [19].However, due to the high computational cost of DEM, these simulations are typically performed on scaled-down geometries and particle quantities [20].The extremely high computational cost of DEM makes this method impractical for design and optimization, because these tasks require several iterations [21].
The amount of highly detailed information yielded by DEM simulations, such as the position and velocity of each particle in the system, is usually unnecessary, for engineers are interested in average properties, such as the mean velocities and volume fractions of the powders.Hence, to extract useful information, users must filter or average the DEM results.These observations suggest that it might be advantageous to formulate transport equations governing the evolution of these mean quantities directly.In this alternative approach, rather than aiming at the detailed solution offered by DEM, one is satisfied with a far reduced description of the flow.Here, each powder is modelled as a continuum that occupies the entire flow domain (all the solid phases interpenetrate each other and the surrounding fluid) [22,23].Therefore, each powder (that is, each solid phase) has an associated set of mass and linear momentum balance equations governing its motion [24][25][26][27][28][29].As for normal fluids, each of these linear momentum balance equations features a (solid) stress tensor, which requires a suitable constitutive equation (or closure equation) able to capture sufficiently well the rheological behaviour of the powder.Since granular flows can exhibit a wide variety of behaviours, in this method (known as Eulerian-Eulerian or multifluid) the key challenge is to develop suitable granular rheological models [30].
While these difficulties do persist, much has been accomplished since Bagnold's pioneering work [31].Various constitutive models have been proposed to describe experimental observations of granular flows [32].In this effort, the first attempts focused on quasi-static deformations and led to elastoplastic yield criteria, such as those of Drucker & Prager [33] and Matsuoka & Nakai (MN) [34], which adopted Mohr-Coulomb failure criteria to predict the failure point for granular materials under load.Notwithstanding, these models do not represent well the flow of granular materials after failure [35].Adopting the Mohr-Coulomb failure principle, Schaeffer [36] developed the first constitutive equation for the stress tensor of granular materials in slow frictional flows.According to his theory, the solid stress tensor is a function of the solid pressure, of the solid rate of deformation tensor and of a constant friction coefficient.This constitutive equation for the solid stress tensor is coupled with the mass and linear momentum balance equations of the solid phase to obtain a mathematical model for the slow frictional flow of the granular material.However, Schaeffer later noticed that the proposed mathematical model is ill-posed and has problematic instabilities [36].Moreover, his model does not consider the effect of particle size on the flow of the powders.A generalized version of his theory relates the value of the friction coefficient μ (that is, the ratio of tangential to normal stresses) to the inertia number I (the ratio of the time scales of the microscopic and macroscopic rearrangements of grains) [37].This more general approach, known as μ(I)-rheology method, was further developed by Jop et al. [38], who extended it into a full tensorial constitutive equation linking the solid deviatoric stress tensor to the rate of strain tensor in the powder via an effective (that is, a non-Newtonian) viscosity.Various methods for increasing the stability and accuracy of the μ(I)-rheology model have been proposed [39,40].So far, the μ(I)rheology model has proven to be the most accurate framework for describing the experimental observations in various granular flows, such as the collapse of granular piles [39,[41][42][43][44].
The present work is part of an ongoing research on the flow and mixing of powders inside continuous horizontal mixers, in particular for pharmaceutical applications [4].To successfully describe the mixing process in these complex systems, a model must first be able to accurately describe the powder flow within them.Developing and testing such a model is indeed a crucial first step towards the reliable simulation of powder mixing in continuous mixers.Accordingly, this paper focuses on the modelling of the flow of monodisperse powders (thus, in this work no mixing is considered).Since the model is intended for industrial applications, the continuum modelling approach is selected.To model the rheology of the granular material, we employ the μ(I)-rheology constitutive equation.As the accuracy of this equation strongly depends on the values of the material properties featuring in it, these values must be evaluated carefully.To this end, we propose a procedure that combines experiments and simulations of the collapse of granular piles.This approach offers a practical, fast, cost-effective, and relatively simple means of determining the required material properties, enabling the application of the μ(I)-rheology model in industrial settings.We also validate the model against (additional) granular collapse experiments.Through qualitative and quantitative comparison of experimental and numerical results, we demonstrate that the μ(I)-rheology constitutive equation is applicable in modelling transient granular flows.The next part of this work is to ensure that the model can successfully describe the flow of monodisperse powders in complex geometries, such as those of continuous horizontal mixers.We start this part by introducing another method for evaluating the material properties of the μ(I)-rheology model, using experimental data of powder flow in a continuous horizontal mixer of industrial interest (Gericke mixer).Then, we validate our numerical results using additional experimental data of powder flow in the same mixer.Typically, the best way to validate a computational fluid dynamics (CFD) model is to compare the field of a variable (for instance, the velocity field) in the experimental and simulation results.But capturing these fields in a continuous horizontal mixer using standard experimental methods, such as PIV, is complex for powders.A far simpler alternative is to compare the experimental and numerical values of a global variable, such as the solid mass loading in the mixer (from here on, we will refer to this variable as mixer loading, denoting it as M).This is a suitable variable for model validation, because it directly relates to the powder rheology.Previous studies have shown that upon a decrease in the feed powder density (i.e. when the bulk density of the powder decreases), the mixer loading increases [4].This counterintuitive behaviour is related to the powder rheology, confirming the suitability of using the mixer loading to validate the results of the simulations.Finally, we use the model to systematically investigate the effect of material properties on the mixer loading.
In performing the simulations of the powder flow in horizontal mixers, we use and compare two numerical techniques: sliding mesh (SM) and multiple reference frame (MRF).While the former is more accurate, it is also considerably more computationally expensive.Thus, having in mind industrial users for the model, we investigate whether the SM method is always needed or the computationally less expensive MRF method is sufficiently accurate.This paper is organized as follows.In Section 2, we discuss the experimental methods.In Section 3, we present the continuum approach and the μ(I)-rheology model.In Section 4, we use the experimental results to validate the CFD model, and we propose a simple method for evaluating the material properties of granular media.Finally, in Section 5, we simulate the flow of powders in a continuous horizontal mixer employing the sliding mesh and multiple reference frame methods, validating the results against experimental data from the literature; moreover, using MRF simulations, we systematically investigate the effect of various material properties on the mixer loading.

Experimental setup and procedure
We carried out a set of collapse experiments of columns of dry powders in a rectangular box made of transparent acrylic plastic plates.To minimize the effects of side walls on the granular flows, we used a box with a width W of 200 mm.In addition, to avoid bending, we used walls with a thickness of 10 mm.Four vertical grooves were excavated on the outer side of the box to fit a vertical metal gate in the channel creating a reservoir for the initial deposition of the powder.The length of the reservoir (set by fitting the metal gate in different locations) could vary between 30 and 100 mm.An illustration of the experimental setup is shown in Fig. S1 of the Supplementary Information.In our experiments, we employed copper powder and quartz sand.Table 1 presents the particle densities, diameters and volume fractions of these powders.For each powder, we defined a set of eight scenarios by varying the aspect ratio of the column a ≡ H i /L i between 0.5 and 6.For each scenario, we filled the reservoir with a mass m i of powder to form a rectangular heap with length L i , height H i and width W of 200 mm.Each reservoir-filling process involved sieving the granular mass and pouring it through a funnel into the reservoir; this resulted in a uniform powder column.We triggered the granular collapse by removing the metal gate suddenly yet carefully.The granular material then spread in the channel until it came to rest.The volume fraction of the powder piles was estimated as ϕ i = m i /ρ s (H i L i W), where ρ s is the solid density of the powder (that is, the density of the solid material out of which the particles are made, which is assumed to be a constant).To ensure the reproducibility of the results, we repeated each experiment three times.

Data acquisition
A digital video camera (Sony α7) was used to capture all the experiments.The experiments were recorded at a rate of 100 frames per second (FPS) to study the dynamic behaviour of the granular flow.The camera was carefully aligned horizontally to ensure that the side view of the collapse was captured (see Fig. S1 in the Supplementary Information).The experimental videos were processed using MATLAB to extract the powder runout.

The fluid dynamic model
This section presents a model for the flow of monodisperse powders containing spherical grains of diameter d and solid density ρ s .We neglect the effect of the interstitial fluid (in our case, air) as it is insignificant [45].The powder flow is described by the velocity vector u(x, t), where x and t denote the position vector and the time, respectively.Assuming that the bulk density of the granular material is a constant (i.e., the granular material is incompressible) and that the particles are packed at the (constant) volume fraction ϕ, the mass and linear momentum balance equations for the granular material read: where g is the gravitational acceleration, ρ ≡ ϕρ s is the powder bulk density, and T is the solid stress tensor.To obtain the velocity and pressure fields we must solve Eqs. ( 1) and ( 2), and for this we require a closure for the tensor field T.

Granular material rheology
In this work, we focus on slow frictional granular flows [37].In such flows, friction between particles is the dominant mechanism generating stress (frictional stress).The slow frictional flow of dense granular materials is similar to that of viscoplastic, generalized Newtonian fluids for two reasons.Firstly, the flow is characterized by a threshold (a yield condition) described by the classical Mohr-Coulomb friction law [22]; this law can be used to determine the failure of the powders and to define a threshold for a switch between the solid-like (when powders resist deformation and maintain their shape under an applied stress; in this state, the powders behave like a solid, with the particles holding their relative positions) and fluid-like (in this state, the powders can easily deform and change their shape, and they flow under the influence of gravity or applied shear forces) behaviours of the powders.Secondly, experimental results have shown that the viscosity of granular materials depends on the shear rate [46].These two observations are taken into account in defining most rheological models for granular materials in slow frictional flow.Note that, in all these models, the slight variations in solid volume fraction that might be present during the flow are neglected; furthermore, it is assumed that the flow of the granular material always remains in the slow frictional regime [38].Drucker and Prager [33] realized that for granular materials the deviatoric stress tensor is pressure-dependent.Based on their yield criterion, Schaeffer [36] proposed a constitutive equation relating the solid stress tensor to the solid (or granular) pressure and the strain rate tensor of the solid phase: where I is the identity tensor, p s is the solid pressure, defined as the mean normal stress (i.e., one-third of the trace of T) and η is an effective viscosity, defined as: where D is the strain rate tensor: and μ is a friction coefficient, a material constant related to the angle of internal friction φ [22]: In Eq. ( 4), |D| is the magnitude (or norm) of the strain rate tensor, defined as: where γ denotes the shear rate.Eq. ( 4) shows that friction increases with the solid pressure and the angle of internal friction.When the shear rate vanishes, the effective viscosity diverges, leading to numerical problems.Moreover, in the simulation of scenarios where the granular medium is initially at rest (for instance, the collapse of granular piles in response to the action of gravity), this rheological model cannot trigger the flowan expected limitation, since this model applies only to flowing powders after they have yielded.To overcome these issues, we introduced a maximum viscosity value, η max .The predictions of the continuum model rely on a thoughtful choice of η max .Small values may lead to unphysical results, while large values may result in numerical instabilities.Here, we used the method proposed by Rauter [47], which was found to agree well with experimental results: √ , H representing a parameter that must be defined for each investigated scenario (for granular collapse flows, this is taken to be the initial height of the granular pile).
The μ(I)-rheology model regards the friction coefficient as a function of the inertia number (also called the normalized strain rate) I.This represents the ratio of two time scales.The first is the microscopic particle relocation time scale (t reloc ≡ d/ ̅̅̅̅̅̅̅̅̅̅ ̅ p s /ρ s √ ), that is, the time scale of particle relocation as a result of applied solid pressure, while the second is the time scale of the deformation of the granular material (t def ≡ 1/γ) [21].Thus, the inertia number is given by: Values of I much lower than unity are found in quasi-static flows, where the macroscopic deformation of the granular medium occurs at a much slower pace compared to the microscopic rearrangement of the grains.Conversely, values of I much greater than unity correspond to flows where the granular medium deforms significantly over a much shorter time than the time required by the grains to rearrange themselves microscopically [46].From simple granular shear flow tests on rough inclined surfaces, Jop et al. [38] concluded that μ(I) takes on a value μ s at vanishingly small values of I and increases asymptotically towards a value μ d at high values of I. Therefore, they proposed: where I 0 is a material constant.The dependence of the friction coefficient μ on the material properties μ s , μ d and I 0 is illustrated in Fig. S2 in the Supplementary Information.The variation of μ with the inertia number I gives far better results than those given by the Shaeffer closure [21], so this is the approach taken in our work.However, the numerical results do depend on the values of the rheological material properties, which are unique to the specific powder at hand; consequently, these values must be evaluated sufficiently well before the fluid dynamic model can be used predictively.This is particularly important in industrial applications, where powders considerably differ from model powders of academic interest (e.g., glass beads); for industrial powders the values of μ s , μ d and I 0 are usually unknown.Finally, it is worth stressing again that the μ(I)-rheology model is suitable solely for slow frictional flows, being unable to accurately model rapid granular flows (such as those observed in gas-fluidized beds).In these flows, other constitutive models, like those obtained via the kinetic theory of granular gases, must be adopted to effectively describe the rheology of the granular medium and accurately predict its fluid dynamics [48].

Multiphase flows
Even if one neglects the presence of the interstitial fluid in which a granular material is dispersed, the granular material (regarded as a continuous solid phase) still generally forms an interface with the surrounding fluid, turning the granular flow problem into a multiphase flow one [49].Therefore, a multiphase flow model is required to capture the interface between the solid phase and the surrounding fluid (in our case, air).In our model, we used the Volume of Fluid (VOF) method [50] to capture the interface between the powder and the surrounding air.Assuming the powder always remains packed (i.e., the bulk density of the granular material is a constant) and the gas phase can be treated as incompressible, the balance equations are: where the subscript m is used for the mixture of powder and air (in reality, a mixture between the two phases is present only close to the interface between the two phases).The model uses a volume fraction α, which is unity in the bulk of the granular phase and zero in the bulk of air, to capture the interface between the two phases.Values between 0 and 1 are present around the interface [50].The evolution of the volume fraction is governed by the following Eq.[50]: which represents a mass balance equation for the granular phase.In Eqs.(10) and (11), the mixture density ρ m is defined as: where ρ e is the density of the surrounding fluid (here air).Note that since α is not a constant, the mixture of the granular phase and air is compressible, even if both constituent phases are assumed to be incompressible.The deviatoric stress tensor τ m is calculated as follows: (14) where μ e and η are the viscosities of the surrounding fluid and of the granular material (given by the granular rheology model), respectively, and D m is the strain rate tensor for the mixture, defined by an equation analogous to Eq. ( 5).The position of the interface is defined to be at any location where α(x, t) = 0.5 [50].

Experimental results
In this section, we first present the granular collapse experimental results and then use these to test and validate the numerical model.Granular column collapse experiments have been a benchmark in studying transient granular flows [30,[51][52][53][54].The column aspect ratio a (that is, the ratio between the initial height H i and length L i of the column) has been identified as the main parameter affecting the flow characteristics and deposit morphology [55].Generally, two collapse regimes are observed in experiments [56].In these regimes, the final runout, defined as L i , scales differently with a : where L i and L f are the initial and final values of the column length L, respectively, c 1 , c 2 and b are fitting coefficients (reported in various experimental and numerical studies), and a * is the transition point between the two regimes [57].Since a study of the relative impact of various parameters on the granular collapse is beyond the scope of this paper, for further details, we refer the readers to the works of Lagree et al. [30] and Lube et al. [56].Here, we briefly present our experimental results and then use these to evaluate the material properties in the μ(I)-rheology model and to validate the numerical model.
We conducted a set of granular collapse experiments with initial aspect ratio a ranging from 0.5 to 6, using copper and quartz powders.Snapshots of the final deposition of the powders for a granular collapse scenario with an initial aspect ratio of a = 3 are shown in Fig. S3 of the Supplementary Information.Fig. 1(a) reports the final runout as a function of a for both powders.As we see, for both powders the value of L * f increases linearly with a for a < 3 and then follows a power law.The second descriptor of the deposit is the final height, that is, the normalized maximum height of the final deposit, defined as H * f ≡ H f /H i , where H f represents the maximum final height of the pile.For short columns (i. e., for a ≤ 1), H * f equals 1, which implies that during the collapse the powder close to the vertical wall of the reservoir does not move.Nonetheless, for taller columns (i.e., for a > 1), for both powders H * f decreases linearly with a (see Fig. 1(b)).

CFD framework for granular collapse simulations
This section aims to reproduce the experimental results using the mathematical model shown in Section 3. As the material properties in the μ(I)-rheology model (μ s , μ d and I 0 ) are unknown for the powders used in our experiments, we treated them as fitting parameters, carrying out several simulations until the experimental and simulation results matched.As we chose a wide channel to carry out our experiments, sidewall effects on the spreading behaviour of the powder were negligible.Therefore, it was reasonable to carry out 2D simulations.The computational domain was a rectangle with dimensions of 50 cm × 35 cm (L × H).
We used the commercial CFD package Fluent 2021 R1 [58].In the VOF method, the primary phase was ambient air, while the secondary phase was the granular material (an incompressible fluid with density ρ ≡ ϕρ s ).We implemented the μ(I)-rheology model for the viscosity of the secondary phase employing user-defined functions (UDFs).The two phases were immiscible, with no mass transfer occurring at the interface between them.
Since the flow was low-speed and incompressible, we employed the pressure-based solver of Fluent, with the PISO (Pressure-Implicit with Splitting of Operators) algorithm used for the pressure-velocity coupling.For spatial gradient discretization, we used the Least Squares Cell-Based algorithm PRESTO!.To discretize pressure and momentum terms, we adopted second-order upwind algorithms, while to discretize the volume fraction, we used the Geo-Reconstruct algorithm.Temporal discretization was first-order accurate and implicit.To compute the flow variables, we used a maximum of 100 iterations for each time step.To reduce the computational time, we used an adaptive time step with a global Courant number of 0.3.Usually, convergence was attained within the iteration limit by setting the acceptable absolute residuals equal to 10 − 5 .
The no-slip and no-penetration conditions were applied for both phases at the bottom wall.Following Nguyen et al. [59], we assumed that the vertical wall was frictionless, employing the free-slip boundary condition (i.e., the shear stress at the wall was set to zero).A pressure outlet with normal backflow boundary condition was applied at the other two boundaries (which are in contact with ambient air).This boundary condition sets a free outflow (i.e., a zero gradient) unless the velocity vectors point into the domain (inflow).If the inflow is predicted, the volume of fluid function switches to zero so that there is only air inflow.
To obtain mesh independence, we performed simulations for a scenario with an initial aspect ratio of 6.In these simulations, the initial length L i and height H i of the powder column were 3 cm and 18 cm, respectively.In conducting the mesh independence analysis, we focused on the final runout (that is, the value of L i attained at the time when the system reached stationary conditions).A Cartesian mesh with a grid size of Δx was used throughout the whole domain, and the mesh size was varied from Δx = L i /20 to Δx = L i /50.The results of the simulations showed that for mesh sizes smaller than Δx = L i /30, the final runout did not change.Thus, in this section, we ran all the simulations with a mesh size of Δx = L i /30.
Note that, at this research stage, the material properties of the powders were still unknown, so we employed the material properties suggested by Lagree et al. [30]; thus, the results obtained in the mesh independence analysis were not valid yet (that is, they could not be compared with the experimental results).

Evaluation of the material properties
To evaluate the material properties, we first selected the collapse scenario with a = 2.We started with an initial guess for the material properties (μ s , μ d and I 0 ) and (adopting the trial-and-error method) carried out several simulations until we obtained the values of the material properties that could reproduce the experimental results.We selected L * f and H * f as criteria for comparison between experimental and numerical results.We regarded the evaluated values of the material properties as acceptable if the difference between the experimental and numerical results was <5% for both L * f and H * f .In addition, we required that in the simulations the evolution of the powder-air interface had to match the experimental results qualitatively.Once we had obtained satisfactory material properties for this scenario, we reproduced the collapse scenarios with an initial aspect ratio of 3 and 4. If the simulation results were unacceptable for any of these cases, we went back to step one, re-evaluating the values of the material properties.We repeated this loop until we found the values of the material properties that could reproduce all three scenarios, both qualitatively and quantitatively.Fig. 2 shows a block diagram for the evaluation procedure of the material properties.
A final word about the initial guess values required at the beginning of the procedure just described.The accuracy of these values affects the number of simulations required to reach convergence; thus, an educated initial guess is preferable.To this end, the user should refer to the literature to check whether values for powders similar to the one of interest are available.If none are available, one can employ the following criteria.As shown in Fig. S2, I 0 marks the transition of the friction coefficient μ(I) from quasi-static flows, in which I≪1, to flows where the mean deformation of the granular material is more rapid than the grain microscopic rearrangements, in which case I≫1.Consequently, one should expect I 0 to have unit order of magnitude.Experimental findings (including ours) confirm this, but show that I 0 often is closer to the lower bound of the unit order of magnitude interval, with values of around 0.3.Therefore, a guess value of about 0.5 would appear to be reasonable.The material properties μ s and μ d represent the limiting values of the friction coefficient in the limits of very small and very large (compared to unity) values of the inertia number.Their values are system-dependent, but again experimental findings seem to suggest that they should have unit order of magnitude [46,60].

Validation of the CFD results
The next step was validating the CFD results obtained using the values of the material properties estimated with the procedure just described.To this end, we simulated the collapse scenarios with initial aspect ratios of 5 and 6 (in these simulations, the values of the material properties were assigned; that is, they were no longer fitted).The simulation results of these scenarios are in the acceptable range for both characteristics lengths (i.e., L * f and H * f ), so we concluded that the values estimated for the material properties were valid, and the results of the CFD model were sufficiently accurate for the system at hand.This evaluation and validation procedure was repeated for both powders.The final values of the material properties are given in Table 2.
To illustrate the accuracy of the model, we present a comparison between the experimental and simulation results (which are carried out using the material properties reported in Tables 1 and 2) for the Cu powder collapse with an initial aspect ratio of 6.The free-fall time scale is generally used to characterize the granular collapse; therefore, we normalized the time coordinate using this time scale.Fig. 3(a) compares the evolutions of the powder-air interface during the collapse obtained both experimentally and numerically.The good agreement between the experimental and numerical results indicates that the continuum model accurately predicts the granular fluid dynamics.Minor disagreements are observed at the early stage of the collapse, and we believe this is caused by the effect of the gate removal, which is neglected in the simulations but can affect the initial spreading behaviour of the powder.
As shown in Fig. 3(a), during the collapse the granular domain can be divided into two regions.The first is the quasi-static region in which the frictional viscosity is very high.Here, as a result, the powder behaves almost like a solid, and no deformations and flow are observed.The

Table 2
Values of powder material properties obtained via fitting.See Fig. 2  second region is the flowing region, where the granular viscosity is lower, and the granular material behaves like a non-Newtonian liquid.
The granular collapse starts with a dense, fast-moving free-fall phase that takes about 1t * .As the collapse continues, the quasi-static region begins to develop from the bottom corner on the right-hand side towards the interface during a spreading phase that lasts until about 2t * .The final stage of the granular collapse is the stopping phase, which takes about 1t * , at whose end the flowing region inside the granular medium has disappeared.Fig. 3(b) illustrates the evolution of the powder runout during the collapse.As one can see, the simulation results agree with the experimental results.
To simulate the powder fluid dynamics in complex geometries (e.g., horizontal mixers) with the continuum approach, one must know the material properties of the powder of interest.The experimental and numerical analysis of the collapse of powder columns is a practical, fast and relatively simple tool for evaluating the material properties of powders.In the next section, we will discuss an alternative method for evaluating the material properties using the experimental results of the flow of the powders in an industrial continuous horizontal mixer.Then, we will analyze the flow of the same powders inside the same mixer, investigating the effect of various parameters affecting the powder flow.

CFD framework for the powder flow in a horizontal mixer
For the second part of this work, we now focus on the flow of a monodisperse powder inside a continuous horizontal mixer; for this, we will employ experimental results from the literature.Ideally, the best approach for validating CFD results is to compare the local values of fluid dynamic variables (such as the velocity and pressure fields) with the respective experimental results.However, since we did not have such detailed experimental results, we had to consider a global variable that is rheology dependent (so that it allows us to assess the accuracy of our CFD simulations) and that is easy to measure experimentally.For powder flows in a mixer, the problem is transient; however, if the operating conditions are kept constant, the flow eventually becomes periodic.Consequently, to validate the model predictions, we can operate in terms of average values instead of instantaneous ones.We selected the (mean) mixer loading M as the global variable for calibrating, testing and validating our simulation results.This is a suitable choice, because the value of M is a function of the mixer geometry, of the powder properties (and therefore of the rheology) and of the mixer operating conditions [61].In addition, measuring the mixer loading dynamics is simple experimentally.
To validate the model predictions, we used the experimental results of Vanarase et al. [62], who investigated the effect of operating conditions on the mixer loading in a continuous commercial horizontal mixer (the Gericke mixer).In some experiments, Vanarase et al. studied the effects of changing the shaft velocity and powder inlet flow rate on the steady-state value of the mixer loading.In these experiments, they examined the mixer loading under three shaft velocities: 39 rpm, 162 rpm and 254 rpm.For each velocity, they provided data on the values of the mixer loading corresponding to powder inlet flow rates of 5, 15, 30 and 45 kg/h.In their experiments, Vanarase et al. fed the mixer with two powders: Avicel PH-200 and APAP.Their properties are presented in Table 3.However, since in the final mixture the volume fraction of Avicel PH-200 was about 97%, in our simulations we assumed that the powder flowing inside the mixer was monodisperse.We reproduced the geometry of the Gericke mixer accurately.The impeller of this mixer consists of 12 blades equally spaced along the axial direction (see Fig. 4).The mixer discharge can be equipped with a weir, a semi-circular disk that permits controlling the mixer loading.Experimental results have shown that the mixer loading changes when the angle between the weir diameter and the vertical axis is altered [62].In all our simulations, the position of the weir was the same as in the corresponding experiments of Vanarase et al.; the chord of the weir was angled 20 degrees from the vertical axis.
In the simulations, the first step was to develop the geometry and mesh for our case study.Generally, the results should be more accurate if the computational geometry is very close to the actual geometry of interest.Complex geometries, such as that of the Gericke horizontal mixer, can be resolved by adopting meshes of high complexity; however, complex geometries have the disadvantage of increasing the computational cost.Thus, in generating the mesh, our first step was to simplify the geometry to decrease the complexity of the mesh.Therefore, in the computational geometry we neglected some details of the mixer design that we deemed to be unimportant (for instance, small volumes between parts, bolts and nuts).Additionally, test simulation results showed that the simulation time decreased by changing the inlet position from the upper surface of the mixer shell to the mixer side (opposite to the outlet).Changing the feed location is deemed acceptable, because it affects only the local variable profiles close to the inlet of the mixer but not the values of global variables, such as the mixer loading, if we assume that the mixer is of a reasonable size.To account for the presence of the weir, we divided the mixer discharge section into two semi-circular subsections (Fig. 4) and applied a wall boundary condition to the weir subsection.
Modelling the powder flow in the mixer is a transient problem.Because the shaft and blades rotate inside the mixer, the simulation domain includes moving solid boundaries (i.e., shaft and blades); this complicates the problem considerably.One can simplify the problem by switching the reference frame from that integral with the walls of the mixer to one where the shaft and the blades are immobile while the walls (and the gravity vector) rotate.This approach does not require a dynamic mesh to account for the rotation of the shaft and blades, and consequently it reduces the computational costs significantly.However, with one reference frame, we still have moving boundaries (the walls).A solution to reduce the complexity of the problem is to define two reference frames, one integral with the shaft and blades (thus, in this frame, these are fixed) and one integral with the walls of the mixer (thus, in this frame, these are fixed).This method is referred to as multiple reference frame (MRF) [58].It divides the computational domain of the mixer into two regions: an inner one integral with the shaft and blades and an outer one integral with the walls.In the former, the momentum balance equation is solved using the non-inertial form (as a result of the change in the reference frame, the Coriolis and centrifugal forces appear in the linear momentum balance equation), while in the latter, the momentum balance equation is solved using the usual inertial form.An interface (which is an imaginary surface) separates the two regions; at the interface, one must apply continuity boundary conditions for the absolute velocity field and for the stress field.
Note that the MRF is an approach that freezes the shaft in a fixed position; that is, it does not account for the change in relative position between the two reference frames.Therefore, this method is approximate, because it freezes the mixer parts relative to each other (in reality, the shaft and blades of the mixer move continuously relative to the walls, so the relative positions of the two reference frames change continuously).This method gives the instantaneous solution for the instant where the shaft and blades and the walls (and inlets and outlets) are at the relative positions considered in the simulation.Thus, in a real periodic problem (like mixers), the MRF method solves the problem for just one specific instant.For a global variable like the mixer loading, the accuracy of this method is acceptable, as we will confirm in the following.
Another method that simulates the flow of powders inside a mixer is the sliding mesh (SM) method [63].The same two reference frames (as the MRF method) are used, one being integral with the shaft (rotating zone) and the other being integral with the mixer walls (stationary zone).In this method, the grid of the rotating zone moves rigidly in the mesh relative to the grid of the stationary zone; hence, the position of the shaft and blades is not fixed with respect to that of the walls.Consequently, this method is more accurate than the MRF method.Unfortunately, the SM method is also far more demanding computationally, so SM simulations are unsuitable for design and optimization, because these require several simulations.In light of this, one of the goals of this section is to compare the simulation results for mixer loading obtained with the MRF and SM methods, to see whether the MRF results are sufficiently accurate.
The first step in mesh generation is to break up the domain into two zones: one stationary and one moving.We divided our computational domain into an inner zone, which contains the shaft and the blades, and an outer zone, which contains all the stationary walls.To generate the mesh, we used the Cartesian method [64] for both the inner and outer zones.An illustration of the geometry and mesh is given in Fig. 4.After a set of mesh independence tests, we set the minimum and maximum grid sizes to 1.6 mm and 2.5 mm, respectively.The total number of cells in the grid was ~3.3 million.
In the MRF and SM methods, we set the rotational velocity of the outer domain to 0 rpm for all the simulations.In each simulation scenario, the rotational velocity of the inner domain was set equal to the impeller velocity.The no-slip boundary condition was applied to all the walls in all simulations.The mass-flow-inlet boundary condition was used at the inlet for the mixture.Assigning the mass flux allows the total pressure to be calculated based on the interior solution.The mass flow rate for the air in all the simulations was set to zero, while the solid phase inlet flow rate equalled the experimental values for each simulation case.At the inlet, the mass flow direction was normal to the inlet surface.The pressure-outlet boundary condition was applied at the outlet surface.Moreover, a backflow condition was applied at the outlet, allowing the air phase to reverse flow into the domain.
The PISO pressure-velocity coupling scheme was applied to couple the velocity and pressure fields [58].The Least Squares Cell Based and PRESTO!discretization schemes were used for gradient and pressure discretization, while the momentum and volume fraction terms were discretized using the second-order upwind and Geo-Reconstruct methods.As this study aimed to predict the final loading of the mixer at steady-state conditions, we did not aim to simulate the whole transient filling period.Consequently, to decrease the computational time, in each simulation (unless otherwise stated), we initialized the system with a specific amount of powder deposited uniformly at the bottom of the mixer.We carried out the SM calculations employing adaptive multiphase-specific time steps, with a global Courant number of 0.3.Usually, we achieved convergence within 120 iterations by adopting a tolerance for all the flow variables equal to 10 − 5 .In the MRF simulations, the maximum number of iterations was set to 10 5 .The SM simulations were carried out in parallel on two Intel Xeon Gold 6248 2.50 GHz (80 cores) processors, while we ran the MRF simulations on one Intel Xeon Gold 6140 2.30 GHz processor core.We chose not to run MRF simulations in parallel since the computational costs are significantly lower than those for SL simulations.

Evaluation of the material properties
In this part of the work, the challenge was to estimate the material properties of the powder used by Vanarase et al. [62] for the μ(I)rheology model.Vanarase et al. did not report the values of these properties.Our solution was to use a part of the experimental results of Vanarase et al. to evaluate the material properties (using the trial-anderror method) and then use the rest of the experimental results to validate the predictions of our simulations.Having selected our approach for evaluating the material properties and validating the model, we then focused on the next challenge: computational cost.The process of mixer filling (even if we start our simulations with an initial nonzero loading) takes about 60 s [65].On average, the time required to simulate 5 s of powder flow in the mixer using the SM method (with parallel simulations on 80 HPC cores) is around 30 h.In Section 4, we discussed the importance of selecting reasonable initial guess values for the material properties.Even with a plausible choice of these guess values, evaluating the material properties requires several simulations, so undertaking this task using only the SM method is unfeasible.Thus, to find reasonable initial values for the material properties, we ran these simulations using the MRF method (which takes considerably less time).Note that, at this stage, we were not confident that MRF simulations could estimate the mixer loading accurately, but starting the SM simulations with values of the material properties estimated via the MRF method significantly decreased the number of trial-and-error SM simulations, and thus the total computational cost.
We selected one experimental operational condition, a rotational velocity of 39 rpm and an inlet flow rate of 5 kg/h.Then, we carried out a few MRF simulations to find an initial value for the material properties.In these simulations, we initiated the process by making an educated initial estimation of the material properties.Subsequently, we conducted MRF simulations to determine the mixer loading.At the end of each simulation, we found the relative error between experimental and numerical results, ( M exp − M sim )/ M exp .The acceptable error range for this part of the research was set at ±12%, a reasonable value, considering that in Vanarase et al. [62] the experimental results are reported without error bars and that we had slightly changed the geometry of the mixer.The simulation time is another contributing factor to considering higher error values in mixer simulations (compared to the granular collapse flows, here the geometry of the system, as well as the flow itself, are much more complex, increasing the time required to complete each CFD simulation).Decreasing the error value would significantly increase the number of simulations required to evaluate the material properties.In our simulations, adopting a trial-and-error method, we changed the values of the material properties of the rheological model.After repeating the MRF simulations five times (each time changing the values of the rheological material properties), the mixer loading settled in the defined ±12%range, so we identified the initial guess values for the material properties.Afterwards, we switched to the more accurate SM simulations.
The next step was to simulate the same scenario using the SM method with the values of the material properties obtained from the MRF method.Since we were not interested in the mixer filling process, to reduce the computational costs, we started our simulations with a mixer with an initial loading of powder deposited uniformly in the mixer.At  the start of the SM simulation, the initial mixer loading was set to 0.18 kg.After a set of trial-and-error simulations, in which we changed the material properties to evaluate their correct values, the mixer loading settled in the acceptable ±12% range.This condition was reached after about 50 s (of the SM simulation time) in each simulation.Even so, we continued the simulations up to 100 s to ensure that the mixer loading remained within the acceptable range.Note that to run the SM simulations for 100 s, we had to run our case in parallel for ~600 h.The next step was to confirm that the initial loading of the mixer did not affect the mixer loading results.Therefore, we simulated the same case with a (random) value of the initial loading higher than the final experimental value of the mixer loading.As illustrated in Fig. 5(a), with an initial loading of 0.44 kg, the mixer loading fell within the acceptable range after about 45 s, indicating that the initial loading of the mixer did not alter the final result.Using the same values for the material properties, we simulated three more cases with 15, 30 and 45 kg/h inlet flow rates (with a shaft velocity still equal to 39 rpm); also in all these cases, as shown in Fig. 5, the mixer loading settled in the acceptable range, again confirming that the material properties had been evaluated correctly.Their values are reported in Table 4.

Validation of CFD results
After having estimated the material properties values, we focused on validating the results of the simulations against the experimental data.To validate the CFD results, we used the values of the material properties evaluated in the previous section and compared the results of mixer loading with the corresponding experimental data.Specifically, we focused on the experimental results obtained at a shaft velocity of 162 rpm.To accomplish this, we conducted four simulations, varying the inlet flow rate in each simulation.To decrease the computational cost, we started the simulations with an initial mixer loading of ~50% of the experimental value.Fig. 5(b) compares the experimental and numerical results of the mean mixer loading.The mean mixer loading refers to the average value of the mixer loading once the flow in the mixer has become periodic.The shaded area identifies the ±12% acceptable range for the simulation results.As we can observe, all the simulation results are within the acceptable range; this means that the fluid dynamic model with the μ(I)-rheology closure and the sliding mesh solution method can reproduce the flow of a monodisperse powder in a horizontal mixer with reasonably good accuracy.
On average, each SM simulation case in Fig. 5 takes ~600 h on 80 HPC cores.This is why we decided to simulate the same cases using the MRF approach and compare the results in terms of mixer loading with those obtained with the SM approach.We repeated all eight SM simulation scenarios using the MRF method, keeping all the parameters the same.Fig. 6 compares the error in estimating the mixer loading between the two techniques.Both methods underestimate the mixer loading, but the MRF method is less accurate, with an error ranging between 14% and 18%.Even so, the MRF method can accurately predict the trend of how the mixer loading changes when the operating conditions change.Moreover, knowledge of the time-varying flow fields (which become periodic after a certain number of shaft rotations) may be unnecessary for industrial applications.On average, MRF simulations take ~29 h on a single core.In light of this, we decided to use the MRF simulations to carry out the sensitivity analysis, because in this analysis we were interested only in trends (and also because this analysis requires several simulations, which renders the SM approach unfeasible).
When selecting between the multiple reference frame and the sliding mesh methods in CFD, one must consider the specific problem at hand and the objectives of the simulations.The MRF method offers computational efficiency and simplicity but may sacrifice accuracy.On the other hand, the SM method provides better accuracy and captures more details, especially in unsteady flows, but requires much more computational resources and careful grid generation.Therefore, in industrial applications, the simulation method must be dictated by the simulation aims.If the user is interested only in trends or in performing a preliminary study of the mixer performance, then the most suitable choice is the MRF method.Conversely, if the user seeks highly accurate results, for instance, to finalize the detailed design of a mixer, then the SM method could be the better option.These factors should guide users in selecting the most appropriate method for a particular CFD simulation.

Sensitivity to the material properties
So far, all the discussed results indicate that the μ(I)-rheology model  can capture the dynamics of complex and transient powder flows.In this section, we used the MRF method to investigate the influence of material properties (μ s , μ d , I 0 , particle size, and powder density) on the mixer loading.We performed a systematic simulation campaign to study the mixer loading sensitivity on each material property.We selected the simulation scenario with a shaft velocity of 39 rpm and a powder inlet flow rate of 5 kg/hr as the reference simulation.For this scenario, the MRF result for the mixer loading is 0.311 kg.A set of six simulations was carried out for each material parameter.The effects of changes on each parameter (normalized by its reference value) on the mixer loading (also normalized by its reference value) are presented in Fig. 7.We see that an In light of the sensitivity analysis, we now focus on elucidating a counterintuitive behaviour observed in some experiments, which showed that the mixer loading increased when the bulk density of the powder formulation fed to the mixer was decreased [4].Our sensitivity analysis shows that the mixer loading decreases following a decrease in powder bulk density when all the other parameters are kept constant.The key to justifying the experimental observations of García-Muñoz et al. [4] is to note that in actual industrial processes switching from one powder formulation to another usually changes many, if not all, of the Fig. 6.Comparison between SL and MRF results for the estimated loading of the mixer when it reaches steady state.material properties.Hence, we investigated how the mixer loading changes when a combination of material properties is changed.The simplest form of this analysis is to change two of the rheological properties at a time (keeping all the others fixed).Thus, we started another MRF simulation campaign to see the effect of changing μ s and ρ in the mixer loading.We chose these two variables since the sensitivity analysis shows that they have the highest impact on the powder viscosity.
The results are illustrated in Fig. 8.As expected, by keeping μ s constant and by increasing only the powder density (or alternatively by fixing the powder density and increasing only μ s ), the mixer loading increases.But when we change these material properties simultaneously, the mixer loading can increase or decrease.For instance, by increasing the powder density by 50% and decreasing μ s by 10%, the mixer loading increases by 25% (yellow arrow in Fig. 8).Nevertheless, a 50% increase in powder density accompanied by a 32% decrease in μ s leads to a ~ 8% reduction in the mixer loading (red arrow in Fig. 8).Fig. 8 clearly shows that the mixer loading depends significantly on both material properties.That being said, it is important to stress that a change in powder formulation would lead to a variation in all the rheological properties of the powder, which in turn would change the effective viscosity of the powder.Thus, the powder flow dynamics inside the horizontal mixer would alter, leading to variation in the mixer loading.Our results indicate that, provided the rheological material properties are evaluated accurately, the continuum method adopted in this research can predict the changes in mixer loading with reasonable accuracy.

Conclusions
In this study, we implemented a continuum model for describing the slow frictional flow of monodisperse powders, using the μ(I)-rheology constitutive equation to model their rheology.Initially, we employed granular collapse experiments to develop a fast, simple and inexpensive method for evaluating the values of the rheological material properties of any powder.After, we investigated the flow of collapsing powders to validate the fluid dynamic and μ(I)-rheology models in simulating transient granular flows.The numerical results matched the experimental data very well, confirming the accuracy of the model.
The primary objective of this research was to simulate the flow of monodisperse powders in continuous horizontal mixers using the continuum approach.We assessed the reliability of our model in describing these flows in an industrial mixer with complex geometry (Gericke mixer).To validate the results of the model, we considered a global variable, the solid mass loading in the mixer, comparing numerical results to corresponding experimental data.Using the sliding mesh approach, we successfully reproduced the powder flow in the mixer with high accuracy.Nevertheless, owing to its considerable computational demand, this method is impractical for design and optimization in industrial applications.Here, we aimed to overcome the challenge posed by high computational cost, making the simulations practical for industrial applications.To this end, we conducted simulations using the multiple reference frame method, significantly reducing the computational time (by approximately 95%) while still obtaining reasonably good predictions of the values of the solid mass loading in the mixer.
Finally, to address a question raised in the literature about the variation of the mixer loading in horizontal mixers following a change in the bulk density of the powder fed to the mixer, we conducted a sensitivity analysis of the material properties of the powders employing our model.Previous studies have indicated that as the bulk density of the feed powder decreases, the mixer loading increases [4].This behaviour surprises, but it is important to note that the mixer loading is influenced not only by the bulk density of the powder but also by the other material properties of the powder.In experiments or industrial processes, changing the bulk density of the powder usually means changing all the material properties of the powder, including the μ(I)-rheology material properties.The variation in mixer loading reflects the changes of all these properties, so it is complex to predict.Our simulations did confirm this, rationalizing the counterintuitive experimental observation and demonstrating the support that our model can offer to gain insight into the complex flow of powders taking place inside continuous mixers.
In conclusion, our work has advanced the understanding and the prediction capabilities for the flow of monodisperse powders in continuous horizontal mixers, offering a practical method for estimating the rheological material properties of granular materials and providing a useable CFD model for describing the flow of these materials inside the mixers.The combination of the μ(I)-rheology model and the MRF simulation method presents a promising framework for the development of accurate multifluid models for the description of the mixing process inside continuous horizontal mixers of industrial interest.

Declaration of Competing Interest
The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Luca Mazzei reports financial support was provided by Eli Lilly and Company.

Fig. 1 .
Fig. 1.(a) Column final runout as a function of initial aspect ratio for copper and quartz powders with different densities and particle sizes.(b) Normalized deposit height versus initial aspect ratio for the same powders.

Fig. 2 .
Fig. 2. Block diagram of material properties evaluation and validation procedure for each powder type (errors are calculated based on the difference between the experimental and simulation results of final runout distance and final column height).

Fig. 3 .
Fig. 3. (a) A comparison between the experimental and simulation results of the temporal evolution of the powder column interface during the collapse.Each figure is overlaid by velocity vectors, and the contours show the shear rate field within the column during the collapse.(b) A comparison between the experimental and simulation results of the evolution of the powder column runout, L * ≡ L − L i /L i , where L is the time-dependent powder column runout during the collapse.

Fig. 4 .
Fig. 4. Left: an illustration of the Gericke mixer geometry divided into two domains.Right: the mesh generated for the geometry.

Fig. 5 .
Fig. 5. (a) Comparison between the experimental and simulation results for the mixer loading.The dashed lines refer to the experimental results, and the solid lines refer to the transient simulation results obtained using the sliding mesh approach.(b) Comparison between the experimental and simulation results for the mean mixer loading.The simulation results are the mean value of M once the evolution of the system (particularly the mixer loading) has become periodic.The shaded area identifies a deviation within ±12% of the experimental data.
increase in μ s , μ d , d and ρ increases the mixer loading, while an increase in I 0 has the opposite effect.Based on the μ(I)-rheology model, this behaviour was expected, since increasing μ s , μ d , 1/I 0 , d and ρ increases the effective viscosity of the powder.We used a linear trendline for each graph to show the dependence of the mixer loading on each material property.

Fig. 7 .
Fig. 7. Sensitivity analysis of (a) μ s , (b) μ d , (c) I 0 , (d) particle size, and (e) powder density on the mixer loading.In all the simulations, the mixer shaft velocity is 39 rpm, and the inlet flow rate is 5kg/hr.In each set of simulations, all the parameters (except the parameter being analyzed) are fixed at the reference value.

Table 1
Properties of the powders used in the experiments.
for the evaluation procedure.

Table 4
Evaluated material properties for Avicel PH200 powder.