Cavitation in Pharmaceutical Manufacturing

Therapeutic proteins are used to successfully treat hemophilia, Crohn’s Disease, diabetes, and cancer. Recent product recalls have occurred because of sub-visible particle formation resulting from the inherent instability of proteins. It has been suggested that particle formation is associated with late stage processing steps of filling, shipping, and delivery. Previous works demonstrated that cavitation might occur in therapeutic vials subjected to agitation or accidentally dropped, but that mitigation can be achieved with fluid property manipulation. The goal of this research was to (1) assess the risk of cavitation under common pharmaceutical manufacturing conditions (i.e., pipe contraction and pumps), (2) establish a simple threshold criterion for when particulate will form, and (3) suggest a series of mitigation techniques based on these thresholds. To accomplish these tasks, computational fluid dynamic simulations for a variety of pipe contraction and fluid properties were performed. The results of this research show that reducing the turbulence in a fluid system will reduce the likelihood of cavitation. Additionally, threshold bounds were created that establish a definitive transition at which cavitation will occur.

Page 2 of 23 human use. In this process, the company is not worried about the quantity of therapeutic protein produced, the priority is the quality of the end product used to gain FDA approval. Once a drug is approved for use in humans and available on the market, companies increase production with the goal of producing more drugs more efficiently. To do this, the pharmaceutical solution is pumped at higher speeds and pressures to increase the efficiency of the manufacturing process. Additionally, the piping, pumping, and associated valves used to transport the solution are altered because of the new specifications for larger sizes, higher flow rates, and higher pressures. The significance of the changes in fluid and flow properties when machinery is changed from small to large scale is not evaluated, therefore the retrofitted system has a high potential to cavitate.

Cavitation
Cavitation is a phenomenon where a liquid flash boils to vapor and then collapses. Although there are some proposed methods to harness hydrodynamic cavitation in ways to increase pharmaceutic manufacturing efficiency13,14, cavitation is usually an undesirable phenomenon. The flash boiling that occurs with cavitation can be caused by hydrodynamic, mechanical, or acoustic forces that lower the local pressure to below the critical pressure. Because this is a local phenomenon, once the external force driving the low pressure is removed, pressure returns rapidly to the bulk value and the vapor bubble rapidly collapses. During the collapse of the vapor bubble to the bulk fluid, high velocities are present and collide at a singularity point creating high pressures and large amounts of localized heat. The subsequent dissociation of water can occur, which generates free radicals. These are localized phenomena making the temperature and pressure gradient extremely high at the site of cavitation15. Although the primary key to understanding when cavitation occurs is well defined, the factors that reduce its onset and likelihood are not widely documented.
Traditionally, it is understood that a fluid will cavitate when the vapor pressure of the fluid is equal or greater than the local low pressure. This is true when no turbulence is present. However, when turbulence is present a fluid will cavitate at a higher pressure called the critical pressure. The vapor pressure is computed by the Clausius Clapeyron equation which is a function of the temperature of the fluid, Equation 116. Several methods have been proposed to calculate the critical pressure17. For this application and under these conditions, the most appropriate equation for the critical pressure is derived by Singhal et al.18,Equation 2, because it utilized the turbulent kinetic energy to account for turbulence in the model. This results in a critical pressure equal to the vapor pressure when turbulence is not present and an increase in the critical pressure when turbulence is present.
(1) 195 Where Pvap is the vapor pressure of the fluid, T is the temperature of the fluid, and A, B, and C are experimentally determined constants provided by Poling et al.16, Pc is the critical pressure,  is density, and k is the turbulent kinetic energy (check literature19 for more information about k).

Modeling Cavitation
Due to the complexity of cavitation and the different ways this phenomenon can occur, modeling cavitation is challenging. Not only does modeling cavitation require a detailed understanding of when a flow cavitates to form a vapor bubble, it also requires understanding when a vapor bubble will collapse back into the bulk fluid. For specific applications, codes can be created and applied on an individual basis, but a generic code that applies across all applications of cavitation is much more challenging. The available codes for specific application include cavitation due to hydrofoils, shock waves, turbulence, and high speed propellers20-27. These models are only valid for their specific application. Singhal et. al. has produced several specific models that are appropriate for this application. Singhal's "Full Cavitation Model" utilizes a robust algorithm with generic cavitation application and allows the cavitation point to be specified directly by the user18. The Full Cavitation Model is currently being utilized by industry and commercially available Computational Fluid Dynamics (CFD) software to model water and oil pumps, inducers, impellers, and fuel injection systems18.

Proteins Denaturation
The ionic and covalent bonds holding the molecules in the protein together are very strong, but, in comparison, the secondary bonds holding the shape and structure of the protein are relatively weaker. These secondary bonds can be broken by either chemical or physical instabilities2. Physical instabilities include denaturation, aggregation, precipitation, and absorption2. When these unique bond structures and shapes of therapeutic protein are broken, the proteins become partially or fully unfolded. These unfolded proteins can group together to form large aggregates. These larger aggregates are regulated by the FDA and are blamed for the adverse effects seen in patients3-10. Recent studies have shown that smaller aggregates, in comparison to previous research, still have the potential to cause adverse effect in patients even though they are not regulated by the FDA28.
Several factors can lead to protein unfolding and aggregation such as amino acid sequence, pH level, temperature, and concentration28. Additionally, these proteins have highly hydrophobic and hydrophilic parts. Because of this, when the protein comes in contact with a gas-liquid interface, it is pushed and pulled to the point where they can change shape and unfold. In terms of protein stability, cavitation results in at least two events that are undesired. The large temperature gradients and additional gas-liquid interfaces produced by cavitation are undesirable because of their potential to cause protein degradation.
Page 4 of 23 research has shown therapeutic protein degradation occurs in the process of shipping proteins to the patient or distributer28,30.

Dimensionless Numbers
Classically, the relative importance of geometry, fluid properties, and flow conditions on hydrodynamic phenomena is not conducted via dimensional quantities. Instead non-dimensional terms are utilized to enable the transfer of the information to a broad range of applications as desired in this case.
Diameter ratio can be used to characterize the geometric configuration of a contracting or expanding flow. It is defined as a ratio of minor diameter to major diameter, Equation 3, where Dr is the diameter ratio, Dminor is the minor diameter at the outlet, and Dmajor is the major diameter at the inlet. (3) The Reynolds number (Re)is the most readily used non-dimensional number in fluid dynamics. It provides the relative ratio of inertial to viscous effects. Although frequently used to monitor the transition from laminar to turbulent flow, it also provides a tool for non-dimensionally describing the dynamic or velocity effects. In Equation 4, ρ is the density of the solutions, ν is the mean velocity of the fluid, D is the characteristic length, and μ is the dynamic viscosity of the fluid.

(4)
Cavitation number is a dimensionless number used to characterize the likelihood of a fluid to cavitate via a ratio of pressure difference between the local flow pressure and vapor pressure of the fluid in relation to the dynamic pressure of the system. The Cavitation number is independent of both the flow geometry and turbulence energy in the flow, see Equation 5, where Ca is the Cavitation number and Plocal is the absolute pressure at a point in the fluid. The cavitation number evaluated in this work was the global minimum cavitation number, and Plocal was determined by taking the minimum pressure over value over the entire model. In the geometry modeled, the local minimum occurred near the edge of the flow constriction.
The pressure drop along the flow direction of a pipe is magnified by the addition of a constriction. The pressure recovery coefficient, Cp, quantifies this pressure loss normalized by the dynamic pressure, Equation 6. Pout is the pressure at the outlet and Pin is the pressure of the inlet. None of the previous non-dimensional quantities explored the significance of the maximum principal stress or turbulence on the potential for a flow to cavitate. As this can lower the critical pressure by 2-fold, its importance should likely not be ignored. Thus, it is important to quantify the Kolmogorov length scale and the energy dissipation rate by developing associated dimensionless parameters.
The Kolmogorov length scale is a length that characterizes the smallest eddies in a flow [31]. Energy is dissipated in a fluid flow through the energy cascade of eddy dissipation. In each eddy, there are two smaller eddies that transfer the energy to a smaller and smaller length scale. The cascade continues until the Kolmogorov length scale is reached which is indicated by no more smaller eddies and at which point the energy is transferred into heat which is dissipated by the fluid. This energy dissipation happens at all scales of fluid flow. The Kolmogorov length scale is determined by the viscous forces over energy dissipation per mass, Equation 7, where η is the Kolmogorov length scale, γ is the kinematic viscosity, and ε is the energy dissipation rate.
The energy dissipation of the fluid was determined by the energy introduced in the system per unit mass. This was calculated using Equation 8 derived from Lengsfeld et al.32 where ṁ is the mass flow rate, and m is the mass in the system.

Model Description:
To provide a controlled model that represents a range of pharmaceutical manufacturing processes, with the potential to cavitate, a constricting square-edge channel was created, as in Figure 1. Parameters that should impact the probability of cavitation in a flow are the fluid properties, the flow properties, and the geometric properties, which will be altered to span conditions found in syringe injection and rapid vial filling. A pipe constriction is a good representation of a commonly occurring feature where cavitation is likely to occur, and this feature is present in piping systems, pumps, and syringes. The fluid, flow, and geometric parameters altered were: major diameter (i.e., inlet), Minor diameter (i.e., outlet), Mass flow rate, Viscosity, Constriction edge type, and Density. Cavitation is a function of the localized low pressure which is altered by fluid velocity, turbulence, and principal shear stress. These factors can be affected by fluid properties such as viscosity, geometric properties such as area constriction, and flow properties such as flow rate and inlet turbulence. The diameter ratio was used to represent changes in the geometric area constriction and the mass flow rate was used to represent changes in the flow properties. The results from the simulation were analyzed using non-dimensional numbers to maximize the application of the findings. The computational model proves to be more stable when a mass flow rate is used to control the inlet of the model, rather than a velocity or pressure inlet. For this reason, the density remains constant, so the mass flow rate is the only factor effecting the flow and velocity of the fluid. The density was set to a value equivalent to that of an aqueous solution.

Fluid Model:
To generate the CFD model ANSYS Fluent33 was used. ANSYS Design Modeler and ANSYS Meshing were used to generate the geometric shapes and mesh for each model. Both of these software packages are capable of being run through text commands and are executed through the DOS prompt, allowing the models in each to be manipulated and controlled through a master MATLAB code. The constricting channels were modeled using a 2D axisymmetric algorithm on a HP Z800 Workstation using up to 7 processors. The models ranged in size from 6,000 to 10,000 cells with a maximum cell skewness of 0.36. The model was run in a steady-state format and implemented with a pressure-based implicit volume of fluid solver. A mass transfer between the liquid and gas phases was implemented using the Singhal et. al.18 full cavitation model. The point at which the fluid underwent a phase change was defined by the liquid critical pressure. This solver utilizes the momentum equation (Equation 9), continuity equation (Equation 10), and the energy equation (Equation 11). A realizable k-ε turbulence model was implemented with standard wall functions. This turbulence model was implemented because it was originally created for modeling internal pipe flows similar to the case exhibited in this fluid model33. The turbulence going into the model is considered negligible with respect to the increase in turbulent energy due to the pipe constriction. The fluid properties and flow characteristics were altered for each model and were determined by the master MATLAB code.
Page 7 of 23 (9) Where ⃗ is the velocity vector, ∇ is the derivative in 3 dimensional space, ̿ is the stress tensor, ⃗ is gravity, and ⃗ is external body force.
Where Sm is a mass source term. (11) Where E is the total fluid energy, keff is the effective thermal conductivity of the fluid, T is the temperature, h is the enthalpy, J is the diffusion flux, ̿ is the effective stress tensor and Sh is a volumetric energy source.

MATLAB Code Setup and Data Generation:
An evaluation of the entire design space encompassed by automated pharmaceutical filling machines was utilized to determine the overall risk of cavitation. The fluid properties, flow characteristics, and geometric shape were altered in a systematic manner so each variable could be analyzed independent of the others. Figure 2 is a flow chart outlining the flow of the CFD solver, the mesh generator, and MATLAB codes. The journal files for the meshing and the CFD solvers were altered by MATLAB to encompass the entire design space and the resulting files were run in the command prompt. The Fluent analysis exported text files with the desired outputs from the analysis and MATLAB read the output variables, calculated the desired values, and plotted the results for analysis. The variables perturbed by the MATLAB code are displayed in Table 1. For each of the columns in Table 1, a square-edge, angle, and fillet channels were evaluated producing 280 data points for each channel type. This resulted in 840 data points for the small-scale channel and an additional 840 data points for the large-scale channels. The 1680 different CFD model were then split into cavitating and non-cavitating cases and analyzed.

Data Analysis Ranges
The range of the variables perturbed by the MATLAB code in Table 1 were chosen to totally encompass the entire design and operating space used by syringe discharge and automated vial and syringe filling. In the manufacturing of theoretic proteins pipe constrictions, needles and syringes can result in diameter ratios from 0.7 to 0.1. This range was chosen because the 0.7 diameter ratio is representative of a minor pipe constriction and 0.1 is representative of a constriction from a syringe to a needle. The mass flow rates were altered from 5 to 200 g/s which represents a slow flow rate used to fill syringes and a fast flow rate used to fill vials34. The viscosity of pharmaceutical solutions changes exponentially as a function of solution concentration due to chain entanglement16. Therefore, viscosity needs to be represented accurately to provide information on how solution concentration can alter hydrodynamic phenomena. Using the exponential curve generated by Giarratano et. al.5 the viscosity of the pharmaceutical solution was determined as a function of concentration. The range for protein concentrations ranged from 0 to 20 mg/mL. This resulted in a viscosity range from 2.89 to 30.168 cP. This range essentially represents a range of fluids behavior from water to motor oil. The density was not altered because it was directly related to the mass flow rate. The inlet was modeled as a mass flow inlet, so altering the density would directly alter the flow rate into the system. Therefore, changing the flow rate and density would have the same effect on the system.

Likelihood of Cavitation
To determine if the operating regions of needle expulsion and automated syringe and vial filling are likely to cavitate, the cavitating and non-cavitating data from the likelihood of cavitation column in Table 1 will be analyzed.

Dimensionless Numbers
The data for the square-edge small scale channel was analyzed to determine trends and how to reduce the likelihood of cavitation. To focus on the factors influencing cavitation the minor diameter was held constant and only the square-edge channel was evaluated. This allowed trends to be identified and characterized using non-dimensional and fundamental fluid mechanics numbers.

Edge Constriction Type
To determine if the type of edge has any impact, the mouth of the constriction was altered to a 45 degree angled wall and a smooth fillet wall, Figure 1. The square-edge, 45 degree and fillet channels were each evaluated at the 280 conditions shown in Table 1 column small scale. The results for each channel were compared to determine what effect the edge constriction type has on cavitation.

Geometric Scalability of Diameter
In each of the above models, the minor diameter was held constant and the major diameter was altered to account for the change in diameter ratio. In this section, an additional 280 cases were evaluated for each channel type. The conditions evaluated are shown in Table 1 column large  Page 10 of 23 scale. These data were compared with the previous results of the small scale and together they were used to create threshold bounds for when a fluid flow is likely to cavitate.

Likelihood of Cavitation:
The 546 different conditions representing the operating regions of needle expulsion and automated syringe and vial filling resulted in 101 cavitating cases. Figure 3 shows a scatter plot with of diameter ratio versus the Reynolds number with the cavitating and non-cavitating data separated. The model is classified as "cavitation" if any amount of vapor is present. Therefore, there is no measure of intensity of cavitation, just if it is present. The Reynolds number was calculated at the exit of the minor diameter for each case. Overlaid on this data are bounds of operation for various manufacturing and delivery conditions. The operating regions for the syringe discharge, syringe filling, and vial filling operations were determined by the operating conditions of a variety of Bosch pharmaceutical manufacturing products34. This figure shows that during needle expulsion, cavitation likely does not occur. However, during syringe and vial filling, cavitation is likely to occur based on the available operating conditions. Both syringe and vial filling run higher risks for cavitation when operating at higher fluid throughputs. This figure clearly demonstrates that cavitation is an issue in pharmaceutical manufacturing and cannot be ignored.

Dimensionless Numbers:
To analyze the factors associated with the likelihood of cavitation, a better correlation between the fluid flow and fluid properties would need to be better understood. Figure 3 demonstrates that at higher Reynolds number the flow is more likely to cavitate but does not characterize cavitation as a function of the fluid flow and fluid properties. Using different dimensionless and fundamental fluids number, the flow can be analyzed to separate the cavitation and non-cavitating cases. Once separated, trends can be undertaken, and mitigations strategies can be gathered to reduce the likelihood of cavitation. A cavitating flow is a function of the lowest pressure in the fluid region. The low pressure in the fluid is a function of the fluid flow and fluid properties, which could make it a leading candidate for separating the cavitating and non-cavitating groups. Figure 4 shows that if the Reynolds number is below 1000, the fluid will not cavitate, but the figure does not give further insight into what happens if the flow is above a Reynolds number of 1000. It also shows that if the low pressure of the fluid is above 20 KPa the flow will not cavitate, and if it is below it is likely to cavitate with only a few cases where cavitation does not occur. This plot begins to lead to a separation of the cavitating and non-cavitating groups but does not explore the key issue of what is causing cavitation. Therefore, another dimensionless or fundamental fluids number is needed. The pressure recovery coefficient is a function of the pressure difference and kinetic energy of the fluid which make it a good candidate for separating the cavitating and non-cavitating regions. Figure 5 shows the pressure recovery coefficient versus Reynolds number with the data separated into groups as a function of diameter ratio, w. The data separates nicely with little overlapping regions. However, Figure 6 shows the pressure recovery coefficient as a function of Reynolds Page 12 of 23 number with the data separated into cavitating and non-cavitating groups, with large areas of overlap. The figures begin to show a separation of cavitating and non-cavitating data, but they were not able to describe a universal threshold bound.  The Cavitation number is a dimensionless quantity that describes how likely a flow is to cavitate. The Cavitation number is a function of the low pressure in the fluid and the kinetic energy. The velocity term in the Cavitation number is with respect to the outlet velocity in the minor diameter. Figure 7 plots the Cavitation number as a function of the Reynolds number with the data separated into cavitating and non-cavitating groups. This figure shows that if the Cavitation number is below ~0.2, cavitation will likely occur for the type of design and conditions tested in this work. Under the conditions considered here, there are many cases where the fluid does not cavitate under a cavitation number of 1. A threshold of cavitation is sensitive to the materials and parameters being considered35. The Kolmogorov length scale represents energy dissipation through viscous effects thus providing detailed information on turbulent strength or intensity. Figure 8 shows the Kolmogorov length scale as a function of the Reynolds number with a contour of diameter ratio. A linear relationship between the logarithm of the Kolmogorov length scale and logarithm of the Reynolds number indicates that the Reynolds number and diameter ratio can be captured in the single robust quantity of the Kolmogorov length scale. Figure 9 shows the Kolmogorov length scale as a function of the Reynolds number with cavitating and non-cavitating groups. Again, there are regions of overlapping cavitating and non-cavitating data, but the relationship between Kolmogorov length scale and Reynolds number still hold strong. With the understanding of the relationship between the Kolmogorov length scale and the Reynolds number, it is worthwhile to plot the Cavitation number versus Kolmogorov length scale. Figure 10 shows the Cavitation number versus Kolmogorov length scale with cavitating and non-cavitating groups with a reasonable separation between the cavitating and non-cavitating groups.

Edge Constriction Types:
To this point, all of the data analyzed have been for a square-edge constriction channel. By changing the type of edge constriction, a smoother transition from the larger diameter to the smaller diameter can be made. This would change the amount of turbulence introduced at the construction resulting in a change of critical pressure and low pressure. Together this will change the likelihood for a fluid to cavitate. In addition to the square edge channel, the 45 degrees angle and fillet edge channel were analyzed. Due to the nature of the edge types, the square-edge channel has the most turbulent energy and the fillet channel has the least. The same 280 conditions, Table 1 column small scale, were evaluated for the 45 degrees angle and the fillet channels. Figure 11 shows the data from the square-edge, 45 degree, and fillet channels plotted as the Cavitation number versus Reynolds number. If a cavitation, transition, and non-cavitating regions are created, the fluid can be generalized by the following: The most revealing fact about the edge constriction is how much effect it has on the likelihood of the flow to cavitate. The mean critical pressure does not change much, but the low pressure of the fluid in the channel has a large effect. For the same 280 conditions the square-edge, 45 degree, and fillet channels resulted in 237, 194, and 8 cavitating cases, respectively. A comparison of the three channels is seen in Table 2.  Geometric Scalability of Diameter: For each of the previous models the minor diameter was held constant and the major diameter was altered to accommodate the change in diameter ratio. It is of interest to know if these threshold bounds can be translated to constricting channels of different scales. A constricting square-edge, 45 degree, and fillet channels were analyzed under the conditions of Table 1 column large scale. Figure 12 shows the Cavitation number versus Kolmogorov length scale for the three edge constriction types with data grouped into cavitating and non-cavitating cases. Additionally, the same threshold bounds created from the previous section are applied. Under these conditions the large-scale data fall into the same bounds as the small-scale data. Figure 13 shows the large scale and small-scale data with all three channel types fitting into the same threshold bounds. These data show the first quantifiable threshold bound for a cavitating flow with respect to the Cavitation number.

Discussion and Recommendations for Pharmaceutical Industry:
The application of these findings with respect to pharmaceutical manufacturing demonstrates that cavitation is a problem in vial and syringe filling. However, when developing new manufacturing systems, steps can be taken to reduce the likelihood of cavitation. The largest and perhaps easiest impact would be to add a fillet to every edge, but this is not always feasible, so at minimum a tapered edge should be utilized. This smoothing of the edge type should be applied in piping systems, pumps, and valves.
Next, much more thought should be applied to reducing the amount of turbulence introduced into the system. Finally, the Cavitation number should be calculated and applied to the threshold boundary conditions to determine if the flow is in a cavitating, transition, or non-cavitating region. If these steps are taken when designing new pharmaceutical manufacturing systems, the likelihood of cavitation can be greatly reduced. Proof of this is evident when comparing the number of cavitating cases in the square-edge channels to that in the fluid channels, 237 compared to 8, respectively. Prefilled syringes may experience protein stability problems. Previous studies evaluated the influence of several parameters on protein stability including the effects of surfactants to reduce protein surface adsorption in syringes10, 36 Page 20 of 23 computational work approached the influence of geometrical parameters on the onset of cavitation in syringes, consequently leading to subvisible particle formation and protein aggregation.

Conclusions
This research shows under current pharmaceutical manufacturing conditions, especially in automated vial and syringe filling operations, cavitation is a real risk that needs to be quantified and mitigated. The most universal method in determining the risk of cavitation is by calculating the Cavitation number and determining which region the flow is operating in. Considering the conditions found in syringe injection and rapid vial filling that were simulated in the present study, if the Cavitation number is larger than 0.2, the fluid will not cavitate. For different conditions not being considered here, it is recommended that modeling of cavitation inception be performed following the guidelines provided in this present study. To increase the cavitation number one can remove sharp edges by rounding, or reduce turbulence by reducing the mass flow rate or increasing the viscosity. These conditions most likely occur at the outlet of pumps, valves, and reductions in cross sectional area (i.e. fill needles).