Due to the increasing computational power and availability of advanced commercial, as well as open-source computational fluid dynamics (CFD) software, performing CFD simulations is becoming more accessible and its use has expanded into many industrial sectors. For the pharmaceutical industry, the topic of mixing is especially of interest, since it is a common unit operation in pharmaceutical processes and empirical and semi-empirical correlations for predicting the hydrodynamic characteristics inside a batch mixing vessel do not consider the complete three-dimensional geometry of the vessel and impeller. This is important in processes such as crystallization, dissolution, compounding, and cellular production in bioreactors, and the obtained information can aid in the selection of the mixer geometry, the choice of mixing speed and for achieving comparable scale-up/down operation of laboratory and industrial mixing vessels. The pharmaceutical industry is not keen on publishing their research and the details of their technology, so the available literature is limited. In this mini review, publications, in which CFD simulations of the mixing process were performed for real pharmaceutical cases are presented, the obtained results are discussed and the software, modeling techniques and methods and their comparisons to measurements are highlighted.
Abbreviations: CFD: Computational Fluid Dynamics; RANS: Reynolds-Averaged Navier–Stokes; API: Active Pharmaceutical Ingredient; SST: Shear-Stress Transport; RST: Reynolds-Stress Model; PIV: Particle Image Velocimetry; DPM: Discrete Phase Modeling
A literature review has shown that computational fluid dynamics (CFD) simulations in the pharmaceutical industry can be divided into the following topics: drug delivery and uptake, respiratory drug delivery, drying, and mixing, which is the topic of this review. CFD simulations have become a valuable tool for process design, optimization, and scale-up . In the pharmaceutical industry, stirred-tank reactors are used due to efficient mixing and the possibility of simple scale-up and appropriate mixing in agitated vessels directly relates to the quality of the product for both drug substance, specifically cell culturing and protein purification, and drug product manufacturing . Simulations offer the possibility to obtain a deeper understanding of the evaluated process and can provide information that cannot be easily measured. Furthermore, CFD simulations consider the complete three-dimensional geometry as opposed to more basic correlations and various operational conditions and vessel/mixer geometries can be tested without the need for laboratory and large-scale operation. The presented CFD modeling papers include simulations made in collaboration with or supported by Roche , Sanofi , Sandoz , Krka d.d. , Pfizer Inc. [1,6], LEO Pharma A/S [7,8], Astra Zeneca Pharmaceutical Development R&D , GlaxoSmithKline [10,11] and one performed with the support of an industrial consortium including BASF, ICI, Malvern Instruments, AstraZeneca, GlaxoSmithKline and Pfizer .
Pharmaceutical Mixing Simulations
Ladner et al.  performed k-ε Reynolds-averaged Navier–Stokes (RANS) simulations with Flow 3D software to obtain shear rates in a stirred vessel, since they cannot be easily measured at specific locations. Shear stress can damage sensitive microorganisms and the cell culture and consequently cause a decrease in productivity and can also influence the active pharmaceutical ingredient (API). Bottom-mounted magnetic stirrers are commonly used in drug product manufacturing, because it is considered that they provide gentle mixing conditions. Magnetic stirrers have a thin gap between the stirrer head and the spigot through which the liquid also flows. They compared the results to experiments with conductivity measurements. They found that although the shear rates in the gap were higher, the high shear rates in the spigot gap can be neglected due to the insignificant flow in that area. They also found that the direction of flow is dependent on the vessel volume, which should carefully be considered when designing small-scale (scale-down) models. Low shear inducing impellers are also used in bioreactors for cell production.
For providing low-shear mixing, the impeller blades usually have a low angle regarding the horizontal position, lower mixing speeds are used, and in the industrial production typically 2 impeller heads can be used for efficient mixing of the larger volume. Ebrahimi et al.  analyzed the effects of this impeller configuration and rotational speed on the mixing performance of a double-impeller bioreactor. As expected, an increase in the mixing speed increased total power consumption and cell stress, while the mixing time was decreased. The CFD model-obtained power numbers were compared to experimentally determined ones from torque measurements on the impeller and a very close agreement was found, which validated their modeling approach, using ANSYS Fluent software and k–ε RANS turbulent equations The typical bioreactor and the CFD solution presenting velocity contours and vectors is presented on Figure 1. In another work, a stirred pilotscale bioreactor was studied by Bach et al.  for Trichoderma reesei fermentations with the standard k–ε RANS model in ANSYS CFX.
The mixing time was determined with tracer experiments in the fermentation broth. The method of using CFD to predict the mass transfer coefficient was validated with experiments and it was also proven that is was as accurate as empirical correlations. Eulerian–Eulerian multiphase flow simulations with the volume-offluid method for capturing the gas-liquid interface were performed by Haque et al.  with the standard k–ε turbulent model, shearstress transport (SST) model and the differential Reynolds-stress model (RST) in ANSYS CFX. The results of the models were similar and comparable to the Laser Doppler velocimetry measurements of the tangential velocity profile. Witz et al.  performed two-phase Euler-Lagrange simulations in large bioreactors with volumes up to 40 m3, which was made possible due to the lattice Boltzmann method applied, which can be highly parallelized and calculated on the many processing cores of modern graphics cards with CUDA technology. Bubble coalescence and breakup were accounted for and a bubble size distribution was obtained, as well as mass transfer coefficients, which were predicted quite well considering the complexity of the simulation, when compared to the measured values (Figure 1).
Pohar et al.  performed simulations of the industrial process of amlodipine maleate cooling crystallization. Apart from heat transfer and population balance modeling they performed k–ε turbulent mixing simulations in Fluent with the rotating frame of reference technique. The concentration of dissolved substance was followed on-line with FTIR measurements as well as the particle size distribution, which was measured with the FBRM probe. They found that the concentration and temperature distribution inside the industrial vessel were homogeneous and good predictions of the particle size distributions under various operation conditions were obtained. The CFD simulation results during cooling crystallization are presented on Figure 2. Another industrial crystallization (combined evaporative/antisolvent) of an API was investigated by studying the effect of hydrodynamics on the agglomeration of particles by Falk et al. . The name of the compound was not revealed but was designated “Compound A”.
Experiments were done at the laboratory scale and it was found that the type of agglomerate formed correlated with the level of agitation. It was presumed that the obtained agglomerate type was dependent on the differences in the collisions of primary crystals due to differences in the local degree of agitation. They found a correlation between the turbulence dissipation rate (ε) of crystals and agglomeration rate into a specific morphology. A flaky, loose agglomerate was produced by conditions with lower ε. Turbulence dissipation rate was therefore used as the scaling parameter for pilot plant and commercial scale crystallization operating conditions. For the modeling, they used the Realizable k–ε turbulence model in Fluent. In another work, pharmaceutical antisolvent crystallization of aspirin from ethanol (solvent) and water (antisolvent) was studied by Öner et al. . CFD simulations were used to obtain the mixing characteristics at different agitator speeds.
The fluid dynamics during the feeding of the antisolvent into the crystallizer were modeled by developing a dynamic compartment model and implementing it in the MATLAB/Simulink environment. The changes of solvent concentration, density, and compartmental volume were taken into account by the model while assuming a dynamic flow between compartments during filling. ANSYS CFX software was used with the standard k–ε turbulence model. A study of a batch crystallizer was also done by Chew et al.  in Fluent, who studied paracetamol crystallization using a realizable k–ε turbulence model and large eddy simulation. They compared the operation of a conventional impeller driven batch crystallizer and an oscillatory baffled batch crystallizer and found that the particles precipitated in the latter case were of significantly higher quality. Simulations have shown that the shear rate in the latter case was much higher and responsible for a higher nucleation rate, which provided particles of a smaller size (Figure 2).
Flow of Particles
Particulate flow and mixing behavior in the blending of dry powders were modeled and simulated by an Eulerian-Eulerian multiphase framework by Nguyen et al.  in Fluent. The equations of dense particulate flow in a mixer with a high shear rate were closed with the application of a frictional stress model and the kinetic theory of granular flow. A scalar transport equation was added as the tracer to capture the transient mixing dynamics of the system. A high-speed camera and particle image velocimetry (PIV) were used to experimentally determine the evolution of the movement of the tracer. The imaging technique was processed by MATLAB image toolbox in order to obtain the local concentration of particles. With this model, the main features of granular flow were captured: the bed height and the dominating direction of flow. They found a typical poor mixing region under the impeller at small velocity gradients, where a considerable solid build-up was noticed.
They used a laminar flow model due to the presumed insignificant effect of turbulence. In another work, Waghmare et al.  developed a scale-up semi-empirical correlation using CFD for the drawdown of floating solids based on the average liquid velocity at the free liquid surface. Fluent was used with the RNG k–ε turbulence model and the rotating reference frame method; sliding mesh was considered but was not chosen due to the required time for the simulations. After the solution of the flow field, the particles were tracked with discrete phase modeling (DPM). Fumed silica was added from the top for the experiments and injected uniformly from the top liquid surface in the simulation. The obtained correlation was also extended to industrial-scale reactors with up to 4 m3 of volume with the CFD simulation. The correlation obtained can significantly decrease the time and effort for scale-up to larger production volumes.
While most pharmaceutical processes involve the mixing of Newtonian solutions, Cortada-Garcia et al.  studied the mixing of a non-Newtonian mixture of glycerol and a gel formed of polyethylene glycol and Carbomer, which are used in the manufacturing of non-aqueous toothpastes. Rheological measurements provided a power-law dependence for the viscosity on the shear rate, which was then used in the CFD simulations. The obtained experimental power curves were in good agreement with the modeling results. The flow was laminar, so a turbulence model was not required. They also compared the results of the rotating reference frame and sliding mesh methods and a 2% difference in the values for the impeller torque was obtained at the best mesh qualities, while the time of simulation increased from 40 min to 36 hours. The simulations were performed with ANSYS Fluent software.
A very limited amount of publications from the pharmaceutical industry on mixing simulations are available. While CFD software is becoming increasingly common in their production and research & development departments, publication of their results and technologies is not favored. Nevertheless, the few publications on the topic cover the most frequently addressed subjects, which are: evaluation of the shear rate in bioreactors, prediction of crystal size distribution and crystal morphology with population balance modeling, determination of mixing efficiency, homogenization time and poor mixing locations, scale-up/down with appropriate calculations of the turbulent kinetic energy, dissipation rate, average fluid velocity, power number, power per volume, among others. Reynolds-averaged Navier–Stokes equations (RANS) are most frequently used for turbulent mixing simulations due to their low computational cost at a proven accuracy, as opposed to direct numerical simulations, large eddy simulations, detached eddy simulations, and other methods.
For mixing simulations, the rotating frame of reference technique (or multiple reference frame) is also commonly applied for similar reasons of low computational cost and proven accuracy, as opposed to the more realistic simulation of the mixer rotation with moving (sliding) mesh methods, which require considerably longer computational times. The lattice Boltzmann method has been shown to be efficient for large-volume two-phase flow simulations, which are harder to simulate by the Navier–Stokes equations. The majority of work was performed in ANSYS Fluent (previously only Fluent) or ANSYS CFX with the standard, Realizable or RNG k–ε models. For modeling particle dissolution and crystallization, population balance modeling is the most frequently applied.
- Falk RF, Marziano I, Kougoulos T, Girard KP (2011) Prediction of Agglomerate Type during Scale-Up of a Batch Crystallization Using Computational Fluid Dynamics Models. Org Process Res Dev 15(6): 1297-1304.
- Ebrahimi M, Tamer M, Villegas RM, Chiappetta A, Ein Mozaffari F (2019) Application of CFD to Analyze the Hydrodynamic Behaviour of a Bioreactor with a Double Impeller. Processes 7(10): 694.
- Ladner T, Odenwald S, Kerls K, Zieres G, Boillon A, et al. (2018) CFD Supported Investigation of Shear Induced by Bottom-Mounted Magnetic Stirrer in Monoclonal Antibody Formulation. Pharm Res 35(11): 215.
- Witz C, Treffer D, Hardiman T, Khinast J (2016) Local gas holdup simulation and validation of industrial-scale aerated bioreactors. Chem Eng Sci 152: 636-648.
- Pohar A, Likozar B (2014) Dissolution, Nucleation, Crystal Growth, Crystal Aggregation, and Particle Breakage of Amlodipine Salts: Modeling Crystallization Kinetics and Thermodynamic Equilibrium, Scale-up, and Optimization. Ind Eng Chem Res 53(26): 10762-10774.
- Waghmare Y, Falk R, Graham L, Koganti V (2011) Drawdown of floating solids in stirred tanks: Scale-up study using CFD modeling. Int J Pharm 418(2): 243-253.
- Öner M, Stocks SM, Abildskov J, Sin G (2019) Scale-up Modeling of a Pharmaceutical Antisolvent Crystallization via a Hybrid Method of Computational Fluid Dynamics and Compartmental Modeling 46: 709-714.
- Bach C, Yang J, Larsson H, Stocks SM, Gernaey KV, et al. (2017) Evaluation of mixing and mass transfer in a stirred pilot scale bioreactor utilizing CFD. Chem Eng Sci 171: 19-26.
- Nguyen D, Rasmuson A, Niklasson Björn I, Thalberg K (2014) CFD simulation of transient particle mixing in a high shear mixer. Powder Technol 258: 324-330.
- Cortada Garcia M, Dore V, Mazzei L, Angeli P (2017) Experimental and CFD studies of power consumption in the agitation of highly viscous shear thinning fluids. Chem Eng Res Des 119: 171-182.
- Chew CM, Ristic RI, Dennehy RD, De Yoreo JJ (2004) Crystallization of Paracetamol under Oscillatory Flow Mixing Conditions. Cryst Growth Des 4(5): 1045-1052.
- Haque JN, Mahmud T, Roberts KJ, Liang JK, White G, et al. (2011) Free-surface turbulent flow induced by a Rushton turbine in an unbaffled dish-bottom stirred tank reactor: LDV measurements and CFD simulations. Can J Chem Eng 89(4): 745-753.